AD-AgOJ*^ 

IIIIIIIIH 


(P 


AD 


Solution  of  Maxwell’s  Equations 
Using  the  Time  Domain  Method  of  Moments 

Final  Technical  Report 
Supplement  -  publications  and  reports 


by 


Dr  Michal  Mrozowski 
May  1994 

United  States  Army 


DTfC 

electe  ^ 

MAY  31 1994  '  " 

B 


EUROPEAN  RESEARCH  OFFICE  OF  THE  U  S.  ARMY 


London,  England 


Contract  number:  DAJA45  -  92  -  C  -  0032 


The  Technical  University  of  Gdansk 
Department  of  Electronics 
80-952  Gdansk 
Poland 


•Approved  for  Public  Relerise:  distribution  unlimited 


94-16117 

/  ■inniii 


94  5  27  08^ 


List  of  publications  and  reports 


1.  M.Mrozowski,  "Eigenfunction  expansion  techniques  in  time  domain”,  Technical  University 
of  Gdansk,  (internal  report),  April  199d. 

2.  M.Mrozowski,  ”  Function  expansion  edgorithms  for  the  time  domain  analysis  of  shielded  struc¬ 
tures  supporting  electromagnetic  waves”,  J.  of  Numerical  Modelling,  April  1994. 

3.  M.  Mrozowski,  ”  Stability  of  the  Time  Domain  Total  Eigenfunction  Expansion  Algorithm”, 
Proc.  1993  Asia-pacific  Microwave  Conf.,  pp.9-21  -  5-24. 

4.  M.  Mrozowski,  "Derivation  of  Stability  Condition  for  the  Time  Domain  Method  of  Moments 
Algorithms  Using  Functional  Analysis  Approach”,  Time  Domain  Modelling  Workshop,  Berlin 
1993. 

5.  M.Mrozowski,  "A  Hybrid  PEEl-FDTD  Algorithm  for  Accelerated  Time  Domain  Analysis  of 
Electromagnetic  Waves”,  IEEE  Microwave  and  Guided  Wave  Letters,  (submitted). 

6.  M.Mrozowski,  "Stability  Condition  for  the  Explicit  Algorithms  of  the  Time  Domain  Analysis 
of  Maxwell’s  equations”,  IEEE  Microwave  and  Guided  Wave  Letters  (submitted) 


Aooesslon  For 

""rtts  qra&i 

DTiC  ".'j: 

n 

Uuu’ij; 

□ 

Jaa  lU’  loD _ 

!  By - 

!  glgt^lbUtl  Qt?/  .  - 


AveliebUity  QodM 

|A?aii  etid/or 
Sp^iosal 


Dial 


EIGENFUNCTION  EXPANSION  TECHNIQUES 
IN  TIME  DOMAIN 
(internal  report)  -  April  1994. 


Michal  Mrozowski 
Department  of  Electronics, 

The  Technical  University  of  Gdansk 
80-952  Gdansk,  POLAND. 


ABSTRACT 

A  new  approach  to  the  analysis  of  M^ocwelI. equations  in  the  time  domain  is  presented.  The  salient 
feature  of  the  method  discussed  is  that  they  are  based  on  the  expansion  of  unknown  functions  into  infinite 
series  of  entire  domain  or  entire  subdomain  basis  functions  and  the  expansion  coefficients  aire  found  by 
means  of  the  method  of  moments  procedure.  Four  new  explicit  algorithms  are  introduced  and  discussed. 
The  stability  condition  is  derived  in  the  general  way  based  on  the  functional  analysis  techniques.  An  fast 
and  efficient  absorbing  boundary  condition  is  proposed.  The  results  of  the  numerical  tests  are  presented 
showing  that  new  algorithms  can  achieve  greater  speed  than  a  traditional  FDTD  technique. 
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1  THEORY 

1.1  Background 

So  fw  we  have  been  considering  the  frequency  domain  formulations  of  the  internal  electrodynamic  eigen¬ 
value  problems.  The  Fourier  transform  introduces  the  equivalence  between  time  and  frequency.  Therefore 
each  problem  involving  one  of  these  quantities  has  two  rdternative  representations.  The  time  domain  for¬ 
mulation  is  often  used  in  computational  physics  where  its  potential  for  solving  parabolic  and  hyperbolic 
differential  equations  arising  in  diffusion,  convection,  conductivity  or  oscillatory  problems  was  recognized 
many  years  ago  [18]  but  in  electromagnetics  the  frequency  domain  formulations  have  traditionally  been 
preferred  for  the  analysis  of  Maxwell  equations.  In  recent  years  time  domain  techniques  have  increasingly 
been  gaining  audience  tdso  in  the  computational  electromagnetics  community.  Time  domain  adgorithms 
solve  initial  value  problems,  but  this  does  not  mean  that  eigenvalue  problems  are  not  covered  by  time 
domain  formulations  because  eigenvalues  can  be  extracted  from  the  time  response  through  signal  pro¬ 
cessing  techniques.  One  advantage  of  the  time  domain  methods  is  that  they  allow  the  investigation  of 
transient  phenomena.  Additionally,  they  allow  in  principle  a  simple  and  wideband  treatment  of  nonlinear 
and  nonstationary  processes.  The  disadvantage  of  the  time  domain  approach  is  associated  mainly  with 
the  difficulties  in  dealing  with  frequency  dispersive  materials. 

The  time  domain  analysis  can  be  carried  out  based  either  on  differential  or  integral  formulations. 
The  integral  formulations,  which  are  particularly  attractive  for  scattering  and  radiation  problems,  suffer 
from  the  fact  that  their  numerical  solution  is  usually  difficult  and  often  atccumulates  discretization  errors 
which  in  turn  leads  to  instabilities  [1].  For  certain  geometries  of  the  scatterers  it  was  possible  to  overcome 
these  problems  and  develop  efficient  algorithms  [2].  However,  the  difficulties  associated  with  the  integral 
formulations  are  absent  in  algorithms  for  solving  differentisd  operators  and  consequently  most  research  in 
the  area  of  time  domain  modelling  of  electromagnetic  fields  was  concerned  with  three  versatile  and  simple 
methods  based  on  differential  representations.  These  three  algorithms  are  known  as  the  Transmission 
Line  Matrix  (TLM),  Spatial  Network  (SN)  and  Finite  Difference  -  Time  Domain  (FDTD)  methods.  All 
three  techniques  are  well  described  in  the  open  literature  [3]-  [9],  books  [32]  and  monographies  [20].  In  the 
TLM  and  SN  methods  the  structure  under  investigation  is  treated  as  a  spatial  network  of  transmission 
lines  and  the  wave  propagation  is  described  by  incident  and  reflected  voltage  impulses  on  the  mesh  lines. 
In  the  standard  FDTD  method  the  Maxwell’s  equations  are  discretized  both  in  space  and  time,  the 
derivatives  are  computed  by  means  of  central  difference  scheme  and  the  fields  are  evaluated  at  discrete 
nodes.  Although  FDTD  and  TLM  (SN)  methods  use  different  concepts  and  are  concerned  with  different 
physical  quantities,  recent  study  by  Celuch  and  Gwarek  [10]  shows  that  each  of  them  can  be  obtained 
from  the  other  by  a  sequence  of  suitable  transformations.  Hence,  time  domain  methods  currently  used 
in  practice  are  formally  equivalent.  Their  salient  feature  is  that  three  dimensional  space  is  discretized 
and  a  mesh  is  formed.  Samples  of  relevant  quantities  at  mesh  points  are  used  to  represent  a  physical 
continuum.  The  sampling  of  the  solution  in  space  is  one  of  possible  forms  of  a  discrete  representation 
of  a  physical  continuum,  characteristic  for  the  approach  which  we  called  domain  division.  Alternative 
technique  is  to  use  series  of  entire  domain  and  entire  subdomain  expansion  functions  to  approximate  the 
field.  Let  us  look  at  time  domain  algorithms  from  the  point  of  view  of  the  alternative  representation  of 
fields  in  order  to  find  formulations  which  may  broaden  the  range  of  options  available  for  the  time  domain 
analysis  of  complex  problems  of  electromagnetics. 

1.2  Series  expansion  with  respect  to  all  space  coordinates 

Let  us  consider  the  structure  shown  in  fig.  1.  It  is  assumed  that  the  inhomogeneity  is  described  by  the 
relative  permittivity  Cr  and  permeability  pr>  both  being  in  general  functions  of  all  three  space  coordinates 
but  frequency  independent.  Under  these  assumptions  the  Maxwell’s  equations  are  given  by: 


P 

P 


- ? - 

«0€r(*,y,2) 

- 

fiOPr{x,y,z) 


(1) 
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Figure  1:  A  cavity  loaded  with  a  dielectric  inhomogeneity 


These  equations  can  be  written  using  the  following  abbreviated  notation 


§9  =  (2) 

where  Li  =  x  (•),  La  =  x  (•)  and  /  and  y  are  vector  functions  representing 

the  electric  and  magnetic  field,  respectively.  The  above  notation  will  be  used  henceforth. 

Replacing  the  time  derivatives  with  their  finite  difference  equivalents  we  get  the  following  marching 
on  time  equations 

/"  =  + 

,"+»/2  =  y"-i/2  +  AtL2/"  (3) 

The  unknown  functions  f,g  are  now  expanded  into  series  of  functions. 

/"  =  53  “"/•(*-  *)  S'"  =  53  y. «)  (4) 

The  set  of  expansion  functions  is  assumed  to  be  complete  amd  the  functions  are  linearly  independent. 
The  functions  are  defined  on  the  entire  domain  or  have  local  supports  but  they  are  in  general  the  time 
independent  functions  of  all  three  spatial  variables.  Substituting  (4)  into  (3)  we  get 

+  (5) 

Taking  the  inner  product  of  (3)  with  the  expansion  functions  /,  and  y,  results  in 

<r,fi>  =  >+A<<Lly"-‘/^/i  > 

<y"'*’‘^*,yt  >  =  <  y"“‘^*,yi  > +Af  <  L2/",y,  >  (6) 

Using  (  5)  the  above  equations  can  be  cast  into  the  following  matrix  form 

o"  =  o"-‘  + 

jn+i/a  _  + 


(7) 
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o 


where  a  and  6  are  the  vectors  containing  expansion  coefficients  and  ^ ,  B,  C,  D  are  matrices  with 
elements  given  by  the  following  inner  products 

^  Bij  =<  L2/J ,  j  > 

Cij  =  <fj.fi>  Dij=<gj,gi>  (8) 

All  matrices  defined  above  are  in  general  dense.  If  the  expansion  functions  are  orthonormal  then  C  and 
D  are  identity  matrices  and  (7)  becomes 

a"  =  a"-*  + 

jn+1/2  _  + Alga"  (9) 

1.3  Series  expansion  with  respect  to  selected  space  coordinates 

So  far  we  have  presented  the  algorithms  in  which  the  time  is  discretized  and  the  function  expansion  is  done 
in  three  dimensions.  We  shall  call  them  total  exprmsion  algorithms.  For  the  case  when  eigenfunctions  of 
a  Laplace  operator  are  used  as  a  basis  we  will  use  the  acronym  TEE  (Total  Eigenfunction  Expansion). 
Another  version  of  the  expansion  algorithm  is  obtained  if  the  discretization  is  in  time  and  one  (or  two) 
selected  coordinate  and  the  expansion  is  done  with  respect  to  two  (or  one)  remaining  coordinates.  We 
shall  present  here  a  case  when  one  space  coordinate  is  discretized.  The  extension  of  this  algorithm  to 
the  case  when  two  space  coordinates  are  discretized  and  the  expansion  is  carried  out  only  for  the  third 
one  is  straightforward.  The  space  is  sliced  into  subdomains  (fig.  2)  and  the  fields  are  expanded  on  each 


Figure  2:  A  structure  discretized  along  one  coordinate 

subdomain  (slice)  into  series  of  expansion  functions.  This  type  of  expansion  we  shall  call  partial  and  for 
the  eigenfunction  expansion  use  the  name  the  Partial  Eigenfunction  Expansion  (PEE).  To  obtain  the  time 
marching  equations  for  such  case  we  have  to  introduce  the  finite  difference  scheme  for  the  calculation  of 
the  space  derivatives  with  respect  to  the  discretized  vairiable.  Suppose  the  structure  was  divided  into  K 
slices  in  the  z  direction  and  the  slices  are  uniformly  spaced  by  the  distance  Ad.  To  enable  the  application 
of  the  central  finite  difference  scheme  for  the  calculation  of  the  derivatives  in  the  z  direction  we  use  the 
technique  used  in  the  FDTD  method.  In  this  technique  [4]  different  field  components  are  defined  on  two 
meshes,  each  off  by  half  the  space  step  from  the  other.  This  arrangement  is  known  as  the  Yee’s  mesh. 
In  the  present  cmatext  Yee’s  mesh  is  given  in  one  spatial  dimension,  which  means  that  relevant  field 
components  are  defined  on  the  slices  that  are  off  by  half  the  space  step.  More  precisely  the  /f  * 

are  given  at  z  =  kAd  and  E*  are  defined  for  z  =  (k  +  l/2)Ad.  Bearing  this  convention  in  mind 
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we  expand  fields  at  a  suitably  situated  slice  according  to 

/?=E  S')  k  =  l...K  (10) 


The  derivatives  in  the  ;  direction  are  approximated  by 


d  ,  -  ipk 

■X— Vt*!/?  »5  7—: 


dz 


Ad 


(11) 


Calculation  of  the  ^-derivatives  involves  operations  on  the  functions  defined  on  the  adjacent  slices.  There¬ 
fore,  instead  of  equations  (2)  we  get 


•^fk  =  Li,gk +  Li,(gk+i  -  9k) 

§^9k  =  LUk  +  LUfk-fk-i) 

Where  we  have  split  operators  Li  and  Lj  into  transverse  and  z-  part  according  to 

Lj.  =  -HtT - (■)  =  kr  ^  ^  (•) 

foc;(*,y)  eoe;(j:,y)Ad 

(•) 


(12) 


tkotkH^.y) 


^  V.  X  (.)  L*,  =  ^ 


Hofii(x,y)^d 


(13) 


In  the  above  equations  by  i  we  have  denoted  a  unit  vector  in  the  z  direction. 

Replacing  time  derivatives  by  the  finite  difference  scheme,  expanding  functions  in  (12)  and  taking 
the  inner  products  for  each  slice  we  arrive  to  the  set  of  equations  similar  to  (7)  with  the  block  diagonal 
coefficient  matrices  given  by 

4  =  qdiagfi:*.^^*]  i  =qdiag[j:*,j;:*] 

£  =  qdiag[£‘]  £=qdiag[£*]  (14) 

The  elements  of  the  submatrices  are  given 

=  <  (I'll  ~  =<  > /ii,  > 

=  <  (I'll  +  I*** )/;■»> Si»  >  =  —  <  I'u/jk-i’ffik  > 

('Ij  =  <  fjt  ■  fit  > 

A*  =  <9jt’9it>  (15) 

For  orthonormal  basis  functions  matrices  Q  and  D_  become  identity  matrices  and  we  get  again  (9). 
Thus  both  total  and  partial  time-domain  expansion  algorithms  are  formally  described  by  the  same  set 
of  equations  (7)  or  (9).  However,  one  important  difference  between  the  total  and  partial  expansion 
algorithms  is  that  in  the  former  the  matrices  (8)  are  in  general  dense  while  in  the  latter  the  matrices  (15) 
are  always  quasi  diagonal. 


1.4  Numerical  cost  of  the  time  domain  expansion  algorithms 

One  drawback  of  the  expansion  algorithms  presented  in  the  previous  section  is  that  they  may  lead  to  higher 
numerical  cost  than  FDTD  and  TLM.  The  inner  products  appearing  in  (8)  and  (15)  are  independent  of 
time  and  can  be  stored  in  look  up  tables.  If  the  cost  of  the  calculation  of  inner  products  is  neglected  then 
the  cost  of  one  time  step  of  the  algorithms  is  determined  by  the  cost  of  matrix  multiplication.  Depending 
on  the  basis  functions  used  (local  support  versus  entire  domain)  the  overaJl  cost  may  vary  considerably. 
Generally  speaking,  assuming  that  expansion  is  done  using  L  functions,  the  cost  of  one  time  step  is 
of  order  O(I’).  If  the  matrices  are  sparse  the  cost  is  lower.  In  the  FDTD  and  TLM  method  with  N 
nodes,  the  numerical  cost  is  of  order  0{N).  Consequently,  expansion  techniques  aie  comparable  in  terms 
of  numerical  cost  and  memory  to  classical  time  domain  algorithms  when  N.  This  condition  will 
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easily  be  fulfilled  in  slightly  and  moderately  perturbed  homogeneous  structures  when  eigenfunctions  of 
an  underlaying  Laplace  operator  are  used  as  the  expansion  basis.  The  expansion  functions  have  then  the 
space  distribution  of  the  relevant  field  components  of  modes  in  the  corresponding  homogeneous  problem. 
Similar  conclusion  can  be  derived  regarding  memory  requirements.  An  important  special  case  is  when  the 
submatrices  are  diagonal.  This  takes  place  for  homogeneous  regions  (or  slices)  when  the  eigenfunction 
expansion  is  used.  The  numerical  cost  is  then  of  order  0(L)  which  means  that  the  expansion  algorithm 
is  faster  than  FDTD  and  TLM  provided  L  <  N.  While  in  the  toted  expansion  technique  this  case  has 
little  practical  importance,  it  can  be  very  effectively  used  in  the  partial  expansion  algorithm  as  shown  in 
Sections  1.5,  and  2.1. 

At  this  point  it  is  important  to  note  that  a  very  efficient  implementation  of  the  expansion  algorithms 
may  be  obtained  if  the  expansion  functions  are  sines  and  cosines.  Let  us  consider  the  toted  expansion 
algorithm  described  by  (9).  Equations(6)  imply  that  at  each  step  one  evaluates  the  inner  products 
<  >  and  <  >  emd  these  inner  products  are  used  to  update  expansion  coefficients. 

For  sine  and  cosine  functions  the  inner  product  for  all  testing  functions  cem  be  computed  in  a  very  efficient 
way  using  the  technique  described  in  [13].  In  this  technique  ^  the  inner  products  are  computed  in  one 
step  in  a  sequence  of  inverse  and  forward  FFTs.  The  procedure  consists  of  three  steps  which  for  the  total 
eigenfunction  expansion  are 

•  using  inverse  3D  FFT  and  and  a"  calculate  /”  and  their  spatial  derivatives 

•  Compute  functions  and  Lj/" 

•  using  forward  3D  FFT  compute  all  inner  products  <  >  and  <  L2f”,gi  > 

Similar  procedure  can  be  used  for  the  partial  eigenfunction  expansion  algorithm  with  the  3D  FFTs 
replaced  with  the  2D  transforms.  The  numerical  cost  of  computations  of  inner  products  using  FFT 
was  discussed  in  detail  in  [13]  and  depends  on  the  number  of  expansion  functions,  number  of  sampling 
points  and  the  inhomogeneity  size  but  generally  it  is  relatively  low  (C>(n log2n),  with  n  proportional 
to  the  number  of  unknowns)  and  additional  advantage  is  that  there  is  no  need  to  create  a  look  up 
matrices  ^  and  ^ .  Again,  if  the  number  of  expansion  functions  is  low  the  time  domain  algorithms  based 
on  the  eigenfunction  expansions  can  be  implemented  using  little  computer  storage  and  with  the  speed 
comparable  to  that  of  the  corresponding  FDTD  or  TLM  methods.  It  should  be  noted  however,  that  for 
large  inhomogeneities,  when  many  basis  functions  are  needed,  the  performance  of  the  function  expansion 
codes  will  be  worse  than  traditional  techniques. 


1.5  Hybrid  expansion  algorithms 

The  discussion  of  the  numerical  cost  involved  in  the  time  domain  method  of  moments  algorithms  implies 
that  the  application  of  the  entire  domain  expansion  functions  makes  sense  only  for  moderate  pertur¬ 
bation.  This  deficiency  can  be  overcome  by  modifying  the  algorithms  so  that  different  techniques  can 
be  used  in  the  different  parts  of  computational  space.  It  seems  that  the  partial  expansion  algorithm 
offers  interesting  possibilities,  especially  with  regard  to  problems  of  scattering  inside  waveguides  (  total 
expansion  algorithms  can  only  be  applied  to  cavities). 

1.5.1  A  hybrid  PEE-FDTD  algorithm 

Consider  combining  the  partial  expansion  method  with  a  classical  FDTD  algorithm.  The  strategy  to 
obtain  an  optimal  algorithm  would  be  following.  Use  FDTD  in  the  regions  in  which  a  fine  resolution 
of  field  is  necessary  (eg.  near  edges,  media  interfaces)  or  large  perturbation  exists  and  apply  partiri 
expansion  in  the  homogeneous  or  slightly  inhomogeneous  regions.  Let  us  consider  the  benefits  which  may 
be  obtained  form  this  approach.  In  principle  any  complete  set  of  functions  can  be  used  in  the  partial 
expansion  algorithms,  but  we  shall  concentrate  on  the  case  when  the  eigenfunction  expansion  (PEE)  is 
used  in  homogeneous  regions.  Let  us  recall  that  in  the  PEE  expansion  functions  are  chosen  in  such  a 
way  that  they  constitute  a  set  of  eigenfunctions  of  the  Laplace  operator  defined  on  2D  region  forming  a 
slice.  In  that  case  each  expansion  function  satisfies  the  boundary  condition  and  field  equations  globally 
over  entire  slice  and  additionally,  because  the  eigenfunctions  of  Laplace  operator  are  orthogonal,  matrices 
A'  * ,  ^  *  and  ^  ^  *  are  diagonal  for  homogeneous  slices  and  the  computations  become  extremely  fast 
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(no  more  thM  9  multiplications  and  18  additions  per  iteration  per  6  expansion  coefficients  for  sine  and 
cosine  basis  functions).  If  the  PEE  were  2dso  used  in  inhomogeneous  slices,  the  speed  of  the  algorithm 
would  decrease.  However  if  a  classical  FDTD  algorithm  is  used  in  inhomogeneity  region  the  numerical 
cost  will  be,  locally,  the  same  as  the  cost  of  an  orthodox  FDTD.  Outside  the  inhomogeneity  the  PEE 
algorithms  ensures  not  only  fast  computation  but  also  minimizes  the  memory  usage.  This  is  because 
the  eigenfunctions  of  the  Laplace  operator  ate  selected  as  the  basis  for  expansion  and  therefore  very  few 
series  terms  at  each  slice  would  be  necessary  to  represent  field  in  homogeneous  region. 

The  concept  of  the  hybrid  algorithm  PEE-FDTD  is  shown  in  fig.3.  The  two  algorithms  are  interfaced 


Figure  3:  The  mesh  arrangement  in  a  hybrid  algorithm 

at  a  common  slice  (say  z  =  z*)  where  the  finite  difference  and  partial  eigenfunction  expansion  meshes 
meet.  The  transition  form  the  FDTD  to  PEE  is  done  in  the  following  way.  Given  a  field  distribution 
at  the  z  =  zt,  provided  by  the  FDTD  part  of  the  algorithm,  the  expansion  coefficients  at  this  slice  are 
found  by  taking  the  inner  product  with  each  basis  functions  of  the  PEE.  The  inner  products  can  be 
calculated  using  numerical  integration  for  a  general  case,  or,  for  Fourier  series  field  representation,  the 
FFT  procedures  cem  be  used.  To  switch  from  PEE  to  FDTD  the  series  (10)  (through  direct  summation 
or  FFT)  arc  calculated  at  the  interfeM:e  plane  at  the  points  required  by  FDTD. 


1.5.2  A  hybrid  scalar  PEE-vector  FDTD  algorithm 

The  time  domain  algorithms  presented  so  far  were  formulated  on  Yee's  mesh  which  is  convenient  for  the 
modeling  of  vector  fields  described  by  Maxwell’s  equation.  The  Maxwell’s  equations  involve  3  electric 
and  3  magnetic  field  components  which  means  that  at  least  6  variables  per  cell  are  needed  to  characterize 
the  fields  at  each  iteration.  Recently,  Aoyagi  and  al  [21]  noted  that  if  the  electric  field  is  divergence  free, 
the  Maxwell’s  equation  can  be  converted  to  a  set  of  three  scalar  second  order  equations.  As  shown  in 
[21]  this  fact  can  be  used  in  a  classical  FDTD  algorithm  to  reduce  both  the  memory  eind  numerical  cost 
of  the  algorithms.  If  a  scalar  wave  approach  is  applied  to  the  PEE  algorithm  even  greater  savings  can 
be  achieved.  For  homogeneous  sections  of  the  structure  the  field  is  a  mixture  of  the  TE  and  TM  modes 
and  as  we  noted  in  the  previous  sections  each  mode  can  be  treated  separately.  The  propagation  of  TE 
and  TM  modes  is  governed  by  a  scalar  wave  equation  which  in  the  time  domain  is  given  by 

+  =  0  (16) 

where  ^  is  a  scalar  potential  whose  transverse  part  satisfies  the  Neumann  or  Dirichlet  boundary  conditions 
and  V  is  the  velocity  of  light  in  the  medium.  It  can  immediately  be  recognized  that  the  transverse  part 
of  a  scalar  potential  cim  be  used  as  the  expansion  term  in  the  PEE  algorithm.  Consequently  for  each 
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eigenfunction  ,  at  each  homogeneous  slice  we  have 

52  52 

where  Sf  is  a  transverse  wavenumber.  The  differentiation  with  respect  to  ;  and  t  is  approximated  by  the 
finite  difference  expressions  yielding 

-  20"  +  0"  .. )  =  0  (18) 

This  leads  to  ^ul  explicit  algorithm 

0"+‘  =  [2  +  -  2i>^(Ad)-2A2/]  0"  +  t;2(Ad)-2A2t(08  +  0^^,  J  -  (19) 

Or  in  a  more  compact  form 

=  Ml  +  +  C-. )  -  ‘  (20) 

where  A  and  B  are  coefficients.  It  is  important  to  note  that  each  time  step  requires  only  two  variables 
per  slice  and  is  realized  with  only  5  floating  point  operations  per  slice.  (NB.  Modern  computers  achieve 
comparable  speed  in  addition  and  multiplication  [34]).  A  vector  PEE  algorithm  requires  6  variables  per 
slice  and  27  floating  point  operations  under  the  same  conditions  (homogeneous  slice).  These  parameters 
have  to  be  compared  with  other  3D  implementations  of  time  domain  techniques  which  require,  under 
the  most  favorable  conditions,  at  least  six  variables  per  cell  and  from  21  [21]  to  30  [34]  floating  point 
operations  per  cell  per  iteration.  Even  if  we  disregard  the  fact  that  the  PEE  uses  only  a  low  number  of 
functions  in  the  homogeneous  regions,  the  number  of  floating  point  operations  required  for  one  iteration 
makes  the  scalar  PEE  three  times  more  efficient  in  terms  of  memory  and  from  4  to  6  times  faster  than 
any  other  published  scheme.  Obviously,  in  this  comparison  one  has  to  bear  in  mind  that  the  scalar  PEE 
is  applicable  only  to  homogeneous  slices.  For  the  remaining  parts  of  the  structure  one  has  to  use  the 
vector  finite-difference  time  domain  algorithm.  To  interface  the  two  algorithms  a  common  slice  has  to 
be  chosen  and  either  the  amplitudes  of  scalar  potentials  have  to  be  calculated  from  the  field  vectors  or 
the  other  way  round.  A  common  quantity  (field  or  potential)  is  calculated  in  the  same  manner  as  for  the 
vector  hybrid  algorithm  discussed  in  the  preceding  section. 

1.6  Stability  of  the  time  domain  expansion  algorithms 

The  time  domain  algorithms  discussed  in  this  study  are  explicit,  which  implies  that  the  time  step  can 
not  exceed  a  certdn  value  or  the  computations  become  unstable.  Presently  we  shall  derive  the  stability 
criterion  for  the  the  total  and  partial  expansion  algorithms.  Before  deriving  the  stability  criterion  let  us 
present  problem  (2)  in  the  form  of  a  hyperbolic  equation.  Taking  the  time  derivative  of  the  first  equation 
and  applying  operator  Li  to  the  second  equation  we  obtain. 

(21) 

Putting  L  =  -LiL2  we  get 

^/  +  L/  =  0  (22) 

Equation  (22)  is  hyperbolic  differential  equation.  Stability  of  the  differencing  schemes  for  parabolic  and 
hyperbolic  equations  has  been  studied  in  [17]  and  [35]  in  context  of  classical  finite  difference  representation 

of  differential  operators.  The  theorems  used  in  that  case  are  also  applicable  to  the  present  method.  To 

investigate  the  stability  of  time  marching  algorithms  for  the  hyperbolic  equations  it  is  useful  to  present 
a  problem  in  a  canonical  form 

A^tR^f  +  Af=0  (23) 

The  time  marching  algorithm  for  the  above  problem  is  stable  if  the  following  conditions  are  fulfilled 

[I7](p.386); 


A  =  A*  >  0,  R  =  R*  >  0 

R-4>0 

4 


(24) 

(25) 
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In  other  words,  for  the  time  marching  algorithm  to  be  stable  it  is  sufficient  that  both  operators  A  and 
R  are  self  adjoint  and  positive  and  additionally  the  operator  R  —  0.25A  is  nonnegative.  The  canonical 
form  (23)  is  obtained  from  (22)  by  simply  writing  it  as 

p^I^/  +  L/  =  0  (26) 

Where  I  is  the  identity  operator. 

Comparing  (26)  with  (23)  we  get  R  =  1/A^t  and  A  =  L.  It  can  readily  be  verified  that  operator  L  is 
symmetric  and  positive.  It  suffices  to  verify  the  condition  (25).  This  condition  is  fulfilled  when 


II 


IIM 

4 


(27) 


or 

At  < 

Thus  the  maximal  time  step  in  explicit  time  domain  algorithms  considered  here  depends  on  the  norm  of 
the  operator  L. 

For  a  self-adjoint  bounded  operator  L  defined  in  the  Hilbert  space  H  the  norm  is  defined  as  [36] 


v/ilLil 


(28) 


||L||=  sup  I  <  Lj/,y  >  I  =  |Am„| 
y€H.II»ll=i 


(29) 


where  Ama,  is  the  largest  eigenvalue  of  L. 

Let  us  estimate  the  norm  of  the  operator  L.  To  this  end  we  consider  a  simpler,  one  dimensional 
second  order  operator  defined  over  an  interval  0  <  i  <  / 


F  = 


-6(x) 


(30) 


where  6(x)  >  0  is  a  time  independent  positive  continuous  function  of  x.  Let  V  denote  the  domain  of 
operator  F  and  assume  that  it  allows  only  square  integrable  functions  satisfying  Diricblet  conditions 
at  both  ends  of  the  interval.  Since  the  space  of  squtue  integrable  functions  is  infinite  dimensional,  the 
operator  F  is  unbounded  and  consequently  its  norm  is  infinite.  It  is  necessary  to  calculate  its  norm 
in  finite  dimensional  space.  This  is  what  happens  in  practice  because  we  always  look  for  approximate 
solution  to  the  problem  in  a  form  of  a  linear  combination  of  finite  number  of  basis  functions.  The  finite 
set  of  basis  function  defines  the  approximate  finite  dimensional  subspace  of  original  domain.  Consider 
the  following  truncated  set  of  basis  functions 


i  <  Nm 


(31) 


The  basis  functions  (31)  span  a  finite  dimensional  space  TfjVw  CV  in  which  the  approximate  solution  is 
sought  for.  Now  it  is  easy  to  find  the  upper  bound  of  the  operator. 


I|F||  =  ||6(x)^(  )||  <  ||6„ax^(  )||  =  ||F„ 


(32) 


Where  bmax  is  the  maximal  absolute  value  of  6(x)  over  the  interval  <  0,/  >.  The  eigenvalues  A,-  of 
operator  Fm  are  given  by 


Xi  = 


and  consequently  the  norm  of  F 


l|F||< 


bmarP 


Similar  derivation  can  be  used  for  F^,  a  finite  difference  analogue  of  operator  F  yielding  [17] 

4 


IIFaII  < 


{Arf)^6n 


(33) 

(34) 

(35) 
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where  Sd  is  the  discretization  step. 

At  this  point  we  can  return  to  operator  ||L||.  Let  us  consider  a  cube  Q  with  the  dimensions  /  x  /  x  / 
We  seek  the  approximate  solution  given  a  finite  number  of  expansion  functions  in  the  form  of  normalized 
products  of  sines  and  cosines 

sin  or  cos  i.kKN^  (3C) 

Let  V  denote  the  domain  of  the  original  problem  given  by  equations  (2).  The  basis  functions  (36)  span 
a  finite  dimensional  space  C  “D.  We  calculate  the  upper  bound  of  the  norm  of  operator  ||L||.  Note 

that  ||L||  <  IILnall-  Where 

—  (co/loCmin/^min)  VxVx(-)  (30 

and 

<min  =  inffr(*,y,2),  Pmin  =  Pr{x ,  V,  :)  x,y,z£Q  (38) 

Using  the  same  procedure  as  above  we  find  the  norm  of  operator  L  in  "Hsm 


1|L||<||L„||<.;Lx 


(39) 


where  I’max  =  icoPoCminPmin)~^^^  is  the  maximum  velocity  for  a  plane  wave  in  the  medium  filling  the 
structure.  Using  the  above  estimation  we  get  the  following  stability  condition 


At  < 


21 

Cmax  t^^/3 


(40) 


If  0  is  a  rectangular  prism  with  the  dimensions  axbxl  and  the  upper  bound  for  i,  /,  m  in  the  trigonometric 
expansion  functions  is  Km,  Lm,Nm  then  the  condition  (40)  becomes 


A<  < 


_ 2 _ 


(41) 


For  partial  eigenfunction  expansion  scheme  using  sine  and  cosine  Km  x  basis  functions  over 
rectangular  a  x  6  slices  with  the  discretization  Ad  in  z  direction,  the  stability  condition  derived  using 
(35)  is 

At  <  -  ^  (42) 

WvH^(^)^  +  (s3)^ 

Obviously  for  other  separable  cylindrical  coordinate  systems  with  the  discretization  Ad  in  z  direction, 
the  stability  condition  will  be  determined  by  the  maximal  value  of  the  transverse  eigenvalue  (separation 
constant)  6^^^  as  the  functional  analysis  will  give 

At  < - ■  ^  (43) 

fmor  Y^mor 


For  the  discretization  of  all  three  coordinates  with  steps  Ax,  Ay,  Ad  we  shall  get  the  well  known 
Courant  condition  [20] 

At  <  - ,  ■  —  _ ^  (44) 

At  this  point  it  is  interesting  to  observe  that  the  derivation  described  above  can  be  used  for  investi¬ 
gation  of  the  stability  of  compact  2-D/FDTD  schemes  considered  by  Cangellaris  [19].  In  fact  a  compact 
2-D/FDTD  algorithm  for  the  analysis  of  waveguides  is  a  special  case  of  one  of  the  partial  expansion 
algorithms  presented  in  Section  1.3  with  the  space  having  two  uniformly  discretized  cartesian  coordinates 
X,  y  and  the  variation  in  the  z  direction  given  by  one  basis  function  of  the  form  exp{—j0z).  For  this  case 
we  get 

At  <  - , - 1 - (45) 

t'ma*^(^)*  +  +  (f 
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This  condition  is  identical  as  the  one  given  in  [33]  and  [19]. 

Finally,  let  us  compare  the  Courant  condition  (44),  valid  for  the  classical  3D  FDTD  and  TLM  schemes 
with  the  condition  for  total  and  partial  eigenfunction  expansion  schemes  using  harmonic  function  expan¬ 
sion.  For  simplicity  let  us  assume  that  a  cubical  region  a  x  a  x  a  is  considered  with  cubical  lattice 
Ax  =  Ay  =  Ad  for  FDTD  and  the  number  of  expansion  function  in  each  direction  is  A’.vr  for  the  method 
of  moments  algorithms.  Let  =  a/Ax.  is  then  the  number  of  discretization  points  along  one 
direction.  The  Courant  condition  can  be  rewritten  as 


A<a  < 


1 

Vmax(j^)y/3 


(46) 


Comparing  the  above  result  with  (40)  we  get  the  ratio  of  the  maximum  allowable  time  steps  for  the 
considered  algorithms 

A<a/  tNm 


A<a  2N^ 

It  is  seen  that  the  method  of  moments  edgorithm  can  operat . 
FDTD  scheme  if 


(47) 

with  greater  a  time  step  than  a  classical 


Nm  <  »  0.638  — 

X  Ax 


(48) 


The  identical  condition  is  obtained  for  the  3D  partial  eigenfunction  expansion  scheme  with  sine  euid  cosine 
basis  functions.  Since  the  idea  of  the  time  domain  expansion  algorithms  consists  in  using  a  fewer  number 
of  unknowns  to  represent  the  fields,  the  condition  (48)  is  likely  to  be  fulfilled  in  many  practical  situations. 


1.7  Absorbing  boundary  conditions 

One  of  the  most  important  issues  in  treating  open  problems  is  a  proper  definition  of  the  boundary  condi¬ 
tions  at  the  places  where  no  physical  boundary  exists.  This  problem  occurs  also  in  the  partial  expansion 
algorithm  (the  total  expansion  can  be  used  only  in  completely  closed  structures)  when  the  geometries 
involving  waveguides  have  to  be  modeled.  The  waveguiding  structure  is  assumed  to  be  uniform  in  the 
propagation  direction  and  guided  modes  propagate  along  the  structure  without  attenuation.  Since  the 
structure  is  infinite,  then  the  application  of  an  algorithm  with  the  space  discretization  in  the  propagation 
direction  would  imply  division  into  an  infinite  number  of  slices.  To  minimize  the  time  of  computations 
and  reduce  the  computer  storage,  the  number  of  slices  has  to  be  kept  as  small  as  possible.  The  way  of 
achieving  this  is  to  consider  only  a  section  of  a  waveguide  and  introduce  the  boundary  conditions  at  its 
both  ends  which  simulate  the  infinitely  long  structure. 

The  problem  of  keeping  the  computational  space  as  small  as  possible,  is  one  of  the  most  fundamental 
in  all  time  domain  methods  using  finite  differencing  scheme  in  the  modeling  of  uniform  infinite  structures 
and  is  known  as  the  problem  of  the  Absorbing  Boundary  Condition  (ABC)  or  Radiation  Boundary 
Condition  (RBC).  The  name  Absorbing  Boundary  Condition  comes  from  the  physical  interpretation 
of  the  conditions  which  model  infinitely  long  structures.  If  the  structure  is  uniform  then  the  forward 
propagating  wave  does  not  give  rise  to  any  backward  traveling  waves.  So  if  a  section  of  the  waveguide 
is  considered  the  conditions  at  its  extremes  have  to  be  specified  in  such  a  way  that  a  wave  impinging 
on  them  can  not  cause  any  reflection.  From  the  point  of  view  of  an  observer  inside  the  guide  the 
wave  would  get  completely  absorbed.  lYaditional  FDTD  schemes  use  usually  different  schemes  based  on 
the  decomposition  of  wave  equation  into  a  operators  describing  one  way  propagation  (towards  or  away 
from  the  boundary),  approximation  of  the  resulting  operator  using  the  Taylor  series  or  Fade  rational 
approximation  and  setting  to  zero  the  expression  representing  the  wave  travelling  back  to  the  interior 
of  the  structure  [20].  This  process  results  in  the  Mur’s  first  and  second  order  conditions  [12].  Mur’s 
conditions  and  their  modifications  are  very  popular  in  time  domain  analysis  of  open  space  problems 
when  the  propagation  is  not  dispersive.  In  the  waveguiding  structures  especially  when  the  dispersion 
is  large  the  Mur’s  ABC  tire  not  satisfactory.  One  approach  used  for  highly  dispersive  structures  is  to 
replace  waveguide  by  a  long  section  of  a  lossy  transmission  line.  Recent  publications  show  [26]  that  with 
this  approach  the  return  loss  better  than  -80dB  can  be  obtained  in  the  octave  band  at  the  numerical 
cost  ranging  from  few  hundred  to  few  thousand  memory  cells  and  floating  point  operations  per  iteration. 
Alternative  approach  sometimes  called  the  time  domain  diakoptics  [22,  23]  is  al»>  being  considered  by 
some  authors. 


Micha/  Mrozowski:  Eigenfunction  expansion  techniques  in  time  domain 


13 


1.7.1  Time  domain  diakoptics 

The  diakoptics  is  the  technique  in  which  a  large  system  is  first  decomposed  into  modules,  each  module  is 
analyzed  separately  and  the  solution  for  the  complex  system  is  achieved  by  finding  the  relation  between 
the  solutions  obtained  for  individual  modules  [29].  Diakoptics,  by  nature  of  the  algorithm,  is  particularly 
well  suited  for  modern  multi  processor,  parallel  and  massively  parallel  computers  [30].  The  concept  of 
system  segmentation  is  frequently  used  in  network  and  field  theory  in  time  and  frequency  domain.  To 
present  this  technique  in  the  context  of  ABC  for  waveguides,  we  shall  express  the  problem  in  the  language 
of  the  theory  of  linear  systems. 

Suppose  the  structure  is  infinite  in  the  z  direction  and  we  want  to  analyze  the  wave  propagation  in  this 
direction  using  the  partial  eigenfunction  expansion  scheme  (PEE).  To  this  end  we  consider  a  section  of 
the  structure  consisting  of  K  slices  spaced  by  Ad.  We  use  Yee’s  mesh  so  that  relevant  fields  are  calculated 
on  the  slices  that  are  off  by  half  the  space  step.  Calculation  of  the  fields  at  last  location  z  =  (K  + 1/2)  Ad 
requires  calculation  of  the  following  expression 

^  ^K+litn)  -  tpKits) 

■^tK+lp{tn)  «  - ^ -  (49) 

This  can  be  done  only  if  we  specify  tl>K+i(in)  which  is  equivalent  to  imposing  an  ABC  at  z  =  (A^  l)Ad. 
To  find  the  vaJue  of  V’A’+i(fn)  consider  a  linear  system  shown  in  fig.4. 


F(8) 


h(tyH(«) 


% 

m 


Figure  4:  A  linear  system 

The  system  has  the  transfer  function  H{s)  where  s  the  variable  of  Laplace  transform.  The  Laplace 
transforms  of  the  input  F{s)  and  output  R{s)  signals  are  related  via 

R(s)  =  Fis)H[s)  (50) 

This  relation  can  be  written  in  the  time  domain  as 

r{t)  =  m*h{t)  (51) 

where  *  is  the  symbol  for  the  convolution  integral  and  h(t)  is  the  impulse  response  of  the  system.  For 
functions  sampled  at  time  instances  lAf,  the  response  r„  =  r(nAt)  is  eveduated  by  the  summation 

F»-l 

r’n  —  At  /n— mhm  (o2) 

m=l 

If  we  assume  that  r(t)  and  f{t)  are  the  fields  at  z  =  ( A+ 1)  Ad  and  z  =  KAd,  respectively  then  using  (51 ) 
we  may  calculate  the  value  of  the  field  at  z  =  (A  +  l)Ad  hence  imposing  an  ABC  at  the  extreme  slice  of 
the  guide.  In  the  case  of  the  PEE  analysis  of  a  waveguide  the  transfer  function  approach  allows  us  to  find 
the  ABC  by  convolving  the  impulse  response  of  an  elementary  slice  of  the  guide  of  the  length  Ad  with 
the  field  at  z  =  Ad.  Referring  back  to  the  concept  of  diakoptics  we  see  that  the  interior  of  the  analyzed 
guide  and  the  ABC  are  two  modules  which  can  be  analyzed  separately  and  the  response  of  the  whole 
system  will  be  given,  in  time  domain,  by  the  convolution.  For  TLM  method  the  time  domain  diakoptics 
treatment  of  ABC  condition  was  initially  implemented  by  individually  convolving  signals  at  each  point 
of  discretization  [23].  This  approach  is  computationally  inefficient  so  very  recently  it  has  been  refined 
for  the  ABC  in  homogeneous  guides  by  first  converting  the  discrete  representation  of  fields  into  modes 
and  convolving  each  mode  separately  [24,  25].  It  has  to  be  noted  that  for  the  PEE  the  eigenfunction 
representation  is  an  inherent  feature  of  the  algorithm  so  the  solution  of  the  ABC  problem  using  time 
domain  diakoptics  seems  to  be  a  good  choice. 
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1.7.2  Low  cost  algorithms  for  efficient  implementation  of  absorbing  boundary  condition  in 

the  PEE 

At  this  point  it  has  to  be  noted  that  unless  the  impulse  response  is  short  the  calculation  of  the  convolution 
is  very  time  and  memory  consuming.  Typically  around  1000  memory  cells  and  a  similar  number  of  floating 
point  operations  at  each  iteration  is  needed  to  obtain  low  return  loss.  In  order  to  reduce  the  computer 
time  it  is  necessary  to  find  algorithms  which  allow  the  computation  of  convolutions  required  to  model 
ABC  fast  enough. 

Finite  impulse  response  approach  The  most  straightforward  approach  to  accelerate  the  computa¬ 
tion  of  the  convolution  is  to  limit  the  duration  of  the  impulse  response.  In  this  way  an  elementary  section 
of  the  guide  is  modelled  via  an  approximate  "black  box”  with  a  finite  impulse  response  FIR.  Bearing 

in  mind  that  the  Laplace  and  Fourier  transforms  are  closely  related  via  the  substitution  s  =  ju,  we  may 

find  the  transfer  function  of  a  section  of  a  waveguide  the  length  Ad  by  finding  its  frequency  response. 
For  the  time  harmonic  (expjui)  fields  the  propagation  in  the  2-direction  is  described  (in  the  frequency 
domain)  by  the  factor  exp[-j0(ij)z]  The  transfer  function  H(ju)  is  then 

H{s)  =  /f(iui)  = 

where  0{ui)  is  the  propagation  constant  and  c  is  the  velocity  of  the  light.  If  the  guide  is  nondispersive 
then  /9(w)  =  Aui/c,  with  A  being  a  constant,  and  hence 

h{t)  =  6(t  -  AAr)  (54) 

where  i(f)  is  the  Dirac  function.  From  the  above  equation  it  can  be  concluded  that  a  nondispersive  guide 
behaves  as  an  ideal,  nondistorting  all  p.'ss  filter  which  introduces  only  the  time  delay  AAr.  If  the  time 
step  At  in  the  PEE  algorithm  is  chosen  so  that  lAt  =  AAr,  with  I  being  an  integer,  then  the  ABC  is 
determined  by  simply  substituting 

^K+l{tn)  =  ^K(tn-l)  (55) 

In  other  words,  in  a  nondispersive  guide  the  field  at  2  =  (A  l)Ad  is  completely  determined  by  the  field 
at  2  =  KAd  by  /  =  AAd/{cAt)  time  instances  earlier.  In  the  dispersive  guide  the  situation  is  different 
because  in  this  case  the  impulse  response  may  be  in  general  infinite.  If  the  frequency  behavior  of  is 
known  then  an  approximate  ABC  ctm  be  obtained  by  taking  the  inverse  Fourier  transform  of  the  transfer 
function  and  time  limiting  the  impulse  response,  so  that  the  convolution  can  be  performed  fast  enough. 

In  order  to  illustrate  this  concept  in  more  detail  let  us  consider  the  electromagnetic  wave  propagation 
in  a  rectangular  waveguide  with  a  local  inhomogeneity.  The  guide  is  excited  by  the  bandwidth  limited 
signrd  with  the  upper  frequency  limit  lower  thrui  the  cutoff  frequency  of  the  first  higher  order  mode.  Under 
these  assumptions  only  a  fundamental  mode  can  propagate  away  from  the  inhomogeneity.  If  we  apply 
the  PEE  method  to  the  analysis  of  this  guide  then  the  exprmsion  coefficients  correspond  to  the  mode 
amplitudes.  The  coupling  between  expansion  functions  (modes)  takes  place  only  at  the  slices  located  in 
the  region  of  inhomogeneity.  Outside  the  inhomogeneity  the  expansion  functions  are  not  coupled  and 
therefore  the  time  evolution  of  the  expansion  function  can  be  found  at  a  marginal  numerical  cost.  Note 
that  since  the  operating  conditions  allow  only  a  single  mode  propagation,  then  if  we  terminate  the  guide 
sufficiently  far  away  from  the  inhomogeneity,  the  higher  order  modes  excited  in  the  inhomogeneity  region 
will  die  out  and  we  will  have  to  specify  only  the  ABC  for  one  propagating  mode  (If  more  than  one  mode 
would  be  allowed  to  propagate  then  a  separate  ABC  would  have  to  be  specified  for  each  propagating 
mode  in  the  same  way  as  for  the  fundamental  mode).  The  dispersion  characteristic  of  modes  supported 
in  the  rectangular  waveguide  are  known  and  consequently  the  impulse  response  function  required  for  the 
ABC  is  also  known 

A(<)  =  (5g) 

where  is  a  symbol  denoting  the  inverse  Fourier  transform  and 

7(w)  =  (57) 

with  e  denoting  the  light  velocity  and  being  the  cutoff  angular  frequency  of  a  mode.  A  cylindrical 
waveguide  is  a  high  pass  filter.  From  the  circuit  theory  we  may  conclude  that  its  impulse  response  is 
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causal.  We  can  now  consider  the  behavior  of  the  impulse  response.  Since  it  is  difficult  to  calculate  the 
impulse  response  analytically  we  replace  the  section  of  the  waveguide  by  an  ideal  nondistorting  high  pass 
filter  with  the  stop  band  (—Ue<^c)-  The  transfer  function  of  an  ideal  filter  is 

//Oa;)  =  (l-G-„.(i^))e->‘‘^^  (58) 

where  Ar  =  Ad/c  and  Cu.(ui)  is  the  gate  function  which  equals  unity  in  the  stop  band  and  vanishes 
outside.  The  impulse  response  can  now  readily  be  found  (normalization  constants  neglected) 

h(t)  =  6(t  -  Ar)  -  (59) 

Note  that  the  second  term’s  roll  off  is  proportional  to  t.  The  distance  between  first  zeros  is  ixjuir. 
and  if  we  truncate  the  impulse  response  at  this  point,  the  ABC  will  be  modelled  by  convolving  using 
approximately  2ir/(u;cAt)  samples.  The  ideal  filter  are  nonphysical  and  therefore  their  impulse  responses 
are  acausal  but  this  analysis  can  be  used  as  a  guideline  for  determining  the  truncation  point  in  real  guides. 

Adaptive  prediction  The  situation  is  more  difficult  when  I3(ui)  is  not  known.  For  this  case  us  write 
the  transfer  function  of  a  slice  of  a  dispersive  guide  as  H(s)  =  Hi(s)Hd{s)  where  Hi{s)  is  nsfer 

function  of  an  nondispersive  guide  and  Hii(s)  represents  the  distortion  introduced  by  the  di  :i.  In 

terms  of  the  "black  boxes”  the  section  of  the  guide  is  now  replaced  by  the  cascade  of  an  all  pa^  londis- 
torting  filter  (a  delay  line)  amd  a  filter  (bandwidth  limited  or  all  pass)  with  nonlinear  phase  characteristic. 
The  Hi(s)  and  are  defined  as  follows 


Figure  5;  A  system  representation  of  a  dispersive  waveguide 


Hi(s)  =  =  (60) 

Hd(s)  =  .ffd(ju/)  =  (61) 

In  a  physical  waveguide  the  dispersion  is  significant  at  low  frequency  and  gradually  vanishes  as  the 
frequency  increases.  Therefore, 

lim  [c)?(w)  —  Aui]  =  0  (62) 

W-^00 

The  value  of  Ar  depends  on  the  discretization  step  but  can  be  arbitrarily  small.  Consequently,  the 
argument  in  the  exponent  may  be  made  small  for  all  w’s  emd  the  transfer  function  HdU^)  can  be 
approximated  by  two  term  Taylor  series 

Hd{s)  =  HiUu)  =  w  1  -  i[c^(u;)  -  AwJAr  (63) 

Let  us  put  P(w)  =  i[cP{u))  —  AvjAr  Since  we  are  considering  only  causal  systems  (/’(w)  satisfies  the 
Paley-Wiener  criterion  [16])  the  P(jui)  is  the  Fourier  transform  of  function  p(<). 

PUi^)=  f  p{t)e~^'^*dt  (64) 

Jo 

If  p(t)  is  discretized  in  the  time  domain  with  the  time  step  At  the  integral  is  replaced  by  the  finite 
complex  Fourier  series 

s 

PU<^) »  (65) 

•=o 


16 


Micb»l  Mrozowski:  Eigenfunction  expansion  techniques  in  time  domain 


The  transfer  function  of  the  system  is  then 

JV 

H{8)  «  HiUu){l  -  PUu))  =  (66) 

isO 

We  can  now  calculate  the  impulse  response  h(t).  The  impulse  response  is  then 

h(t)  =  6(t  -  AAt)  -  ^  piS{t  -  (iAt  +  AAt))  (67) 

i=0  N 

The  impulse  response  given  above  describes  the  cascade  of  a  delay  line  and  a  digital  hlter  of  the  order 
N  with  unknown  weighting  coefficients  pi.  The  unknown  coefficients  can  be  found  using  the  system 
identification  approach.  One  possible  solution  could  be  to  select  two  adjacent  slices  located  at  z  =  kAd 
and  z  —  {k  +  l)Ad,k+  I  <  K  which  have  the  same  cross  section  as  slices  terminating  the  guide.  The 
relation  between  input  and  output  signals  between  these  two  adjacent  slices  is  the  same  as  for  the  A'-th 
and  K  +  1-st  slice.  We  may  now  sample  the  signals  at  the  input  and  output  of  the  test  section  and 
find  the  set  of  coefficients  pi  which  minimizes  the  error  between  the  actual  response  and  the  response 
predicted  by  (67).  An  advantage  of  using  a  digital  6lter  approach  is  the  possibility  to  synchronize  the 
filter  with  the  time  step  of  the  field  modelling  method.  This  allows  one  to  take  into  account  the  effect 
of  numerical  dispersion.  Additionally,  the  coefficients  p,-  can  be  computed  at  a  low  cost  using  recursive, 
self  adapting  algorithms  [15].  It  should  be  mentioned  that  system  identification  models  have  been  very 
recently  used  with  success  to  accelerate  the  computations  of  the  time  domain  electromagnetic  simulators 
[28]  or  diakoptics  in  which  modules  were  modeled  with  classical  FDTD  algorithm  [27]. 

Off  line  identification  The  evaluation  of  the  convolution  can  be  immensely  accelerated  if  the  response 
function  can  be  represented  suitably  chosen  discrete  model.  Such  an  approach  has  been  already  proved 
to  be  very  useful  in  development  of  time  domain  algorithms  for  frequency  dependent  materials  [20]  and 
therefore  it  is  useful  to  consider  this  technique  for  the  realization  of  the  ABC  in  the  PEE  method. 
Suppose  the  impulse  response  is  given  by 

L 

ft(m)  =  '^a,g,(m)  (68) 

1=1 

with  L  and  gi  being  a  small  integer  and  a  function  defined  for  0  <  (  <  oo,  respectively.  Let  us  assume 
that  the  expansion  functions  are  given  by 

g,(m)  =  e-“'’"  (69) 

We  additionally  assume  that  the  expansion  functions  can  be  evaluated  recursively  so  that 

gi{m  -hi)  =  Aigi(m)  (70) 

with  At  =  exp(— ai)  The  convolution  is  expressed  then  as 

n-"l  L 

r(n)  =  Af  ^  aigt{m)f{n  -  m)  (71) 

m=l  1=1 

Using  the  recursion  (70),  the  convolution  is  evaluated  through  a  compact  formula 

L 

r(n)  =  Af53P'(»»)  (72) 

1=1 

where 

p/(n)  =  atAifin  -  1)  +  A/pi(i»  -  1)  (73) 

It  is  seen  that  the  computations  are  reduced  to  merely  an  update  formula  for  pi.  For  the  case  of  frequency 
dependent  material  considered  in  [20],  the  impulse  response  is  specified  analytically  so  the  expansion 
coefficients  could  be  csJculated  using  the  Prony  method.  Alternatively  a  suitable  recursion  algorithm  can 
be  implemented  using  the  off  line  identification  process  similar  to  the  on  line  technique  described  in  the 
preceding  subsection.  The  advantage  of  the  off  line  identification  is  that  the  ABC  can  be  modelled  prior 
to  time  domain  computations  and  hence  reduce  the  numerical  complexity  of  the  PEE. 
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Figure  6:  Geometry  of  a  rectangular  resonator  loaded  with  a  dielectric  slab  used  in  numerical  test  to 
validate  TEE  and  PEE  algorithms.  The  dimensions  are  width  a  =  12Ad,  heigth  6  =  6Ad,  length  c  =  8A/. 
slab  parameter  tr  =  3.75,  slab  centered,  slab  with  w  varykug. 

2  VALIDATION  OF  THE  ALGORITHMS 

2.1  Time  domain  algorithms 

All  time  domain  algorithms  presented  above,  except  for  the  hybrid  scalar  PEE-FDTD  scheme,  have  been 
implemented  and  tested  for  the  case  of  the  rectangular  resonators  or  waveguides  containing  dielectric 
slabs.  The  eigenfunctions  of  the  Laplace  operator  with  Dirichlet  or  Neumann  boundary  conditions  were 
used  as  the  expansion  functions  and  the  inner  products  were  cdculated  using  the  Fast  Fourier  Transfrom 
approach.  The  algorithms  have  been  found  stable  provided  the  time  step  did  not  exceed  the  critical 
values  given  in  Section  1.6.  Some  numerical  results  for  the  expansion  algorithms  with  sine  and  cosine 
expansion  functions  are  given  in  Table  1  where  the  normalized  resonant  frequencies  Ad/A  are  compared 
for  fundamental  mode  in  a  rectamgular  resonator  obtained  from  a  section  of  a  waveguide  shown  in  fig.6 
using  the  total  and  partial  eigenfunction  expansions  algorithms  with  the  data  FD-TD  published  in  [9] 
and  the  transverse  resonance  method.  The  parameter  Ad  is  the  discretization  step  used  in  the  FDTD 
and  PEE  algorithms.  The  agreement  is  very  good  for  2dl  cases  included  in  the  table.  Subsequently  a 


Table  1:  Comparison  of  the  results  for  the  normalized  resonant  frequencies  Ad/ A  obtained  using  expansion 
algorithms  with  the  published  data  [9]  for  the  rectangular  resonator  with  a  dielectric  slab  (fig.6) 


TRM 

Partied 

Expansion 

Toted 

Expemsion 

FD-TD[9] 

ti;  =  0 

0.07511  (exact) 

0.07502 

0.07511 

0.0750 

w  =  2Ad 

0.0522 

0.05253 

0.05153 

0.0517 

w  =  4Ad 

0.0445 

0.04456 

0.04433 

0.0442 

hybrid  vector  PEE-FDTD  algorithm  presented  in  Section  1.5.1  was  tested.  For  this  algorithm  the  cutoff 
frequency  of  the  Elfn  of  a  rectangular  guide  loetded  with  a  dielectric  slab  shown  in  fig.7  was  computed 
and  compared  with  the  results  obtained  with  a  classical  FDTD  technique.  The  dimensions  of  the  guide 
were  taken  to  be  20  by  6  mm.  The  slab  is  placed  on  the  wider  wall  in  the  symmetry  plane  of  the  structure 
and  has  the  width  4mm  and  relative  permittivity  of  2.5.  Three  cases  were  considered  corresponding  to 
the  slab  height  of  0, 4  and  6  mm.  Although  the  structure  is  symmetric,  the  symmetry  was  not  exploited 
during  the  computation.  In  all  tests  the  identical  excitation,  space  discretization  (Ad  =  .5mm,  time  step. 
At  =  .9Ad/(coVS),  and  number  of  samples  AT  =  8000  were  assumed.  The  results  are  ipven  in  Tables  2 
and  3  together  with  the  CPU  time  for  a  DELL  466/L  personal  computer.  For  the  FDTD  algorithm  alone 
the  CPU  time  for  the  assumed  discretization  mesh  was  53s  and  was  identical  for  all  three  cases. 
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Figure  7:  Geometry  of  a  rectangular  guide  loaded  with  a  dielectric  slab  used  to  validate  the  vector  hybrid 
time  domain  method  of  moments  showing  PEE  and  FDTD  meshes 


Table  2;  Comparison  of  the  results  and  CPU  times  for  the  cutoff  frequencies  of  Elf  n  mode  in  the  20  by 
6  mm  rectangular  guide  loaded  with  a  dielectric  slab  (cr  =  2.5,  w  =  4mm,  h  -  varying)  using  only  FDTD 
and  PEE  algorithms  for  number  of  expansion  functions  N  and  the  localization  of  the  interface  planes  as 
parameters 


Slab  height 

FDTD 

CPU 

PEE 

CPU 

5  =  0  mm 

26.0475  GHz 

53s 

26.0625  GHz 

2s 

h  =  6  mm 

19.9075  GHz 

53s 

19.0375  GHz 

2s 

h  =  4  mm 

20.1875  GHz 

53s 

For  the  slab  height  of  0mm  and  6mm  and  the  discretization  along  the  wider  side  of  the  waveguide, 
it  is  possible  to  use  PEE  algorithm  with  homogeneous  slices  only.  Since  only  one  expansion  function 
is  needed  for  the  mode  considered,  the  CPU  time  is  2s  or  26.5  times  faster  than  the  FDTD.  For  the 
slab  height  of  4mm  slices  in  the  slab  region  are  inhomogeneous  and  therefore  the  hybrid  vector  PEE- 
FDTD  algorithm  was  used.  The  slab  region  was  treated  with  the  FDTD  and  the  lateral  homogeneous 
regions  were  calculated  with  PEE.  The  CPU  time  depends  on  the  localization  of  the  interface  plane 
and  the  number  of  expansion  functions  used  in  the  PEE  part  of  the  algorithm.  When  the  interface  is 
at  X]  =  7.5  and  xj  =  12.5mm,  i.e.  when  the  FDTD  mesh  is  terminated  only  one  slice  away  from  the 
inhomogeneous  region,  the  CPU  time  of  the  hybrid  algorithm  varies  from  18s  (for  1  expansion  function) 
to  278  (for  5  expansion  functions),  with  most  of  the  time  (14s)  consumed  by  the  FDTD  computations. 
The  error  introduced  by  low  number  of  expansion  functions  is  the  largest  if  only  one  term  is  used  but 
for  this  structure  is  lower  then  0.1%  compared  the  result  obtained  from  pure  FDTD  calculations.  The 
error  decreases  as  the  interface  plane  is  moved  away  from  the  inhomogeneity.  This  is  due  to  the  fact 
that  higher  order  terms  in  the  field  expansion  correspond  to  higher  order  waves  traveling  in  the  lateral 
direction.  These  modes  are  attenuated  so  their  contribution  to  the  field  of  the  EHn  mode  becomes  less 
significant.  If  the  interface  planes  are  located  at  xi  =  5  and  xj  =  15mm,  only  one  term  in  the  PEE  part 
is  sufficient  to  obtain  exactly  the  same  results  as  with  purely  FDTD  technique.  The  CPU  time  for  this 
case  is  31s  of  which  28s  is  spent  in  the  FDTD  part.  The  increase  in  the  number  of  expansion  terms  does 
not  change  the  results. 

No  tests  so  far  have  been  carried  out  for  the  hybrid  scalar  PEE-FDTD  discussed  in  Section  1.5.2  but 
once  scalar  verrion  of  PEE  it  is  about  5  times  faster  than  a  vector  algorithm  even  larger  speed  up  factors 
may  be  expected  for  the  structure  considered  above. 
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Table  3;  Comparison  of  the  results  and  CPU  times  for  the  cutoff  frequencies  of  EHn  mode  in  the  20 
by  6  mm  rectangular  guide  loaded  with  a  dielectric  slab  fr  =  2.5,  u>  =  4mm,  h  =  4mm  using  a  hybrid 
algorithm  with  the  number  of  expansion  functions  N  and  localization  of  the  interface  planes  as  parameters 


rr 

mmsmmi 

Error 

rel.  to  FDTD 

CPU 

FDTD  part 

CPU 

PEE  part 

CPU 

combined 

Speed  up  rel. 
to  FDTD  (53s) 

Interface  of  algorithms  at  r  =  7.5,  x  =  12.5  mm 

ai 

20.2075 

148 

48 

188 

2.9 

m 

20.1975 

4-0.05% 

14s 

98 

238 

2.3 

Bl 

20.1975 

148 

138 

278 

1.96 

Interface  of  algorithnns  at  x  =  6.5,  x  =  13.5  mm 

HI 

20.195 

19s 

4s 

23s 

2.3 

m 

20.1925 

-1-0.01% 

198 

88 

278 

1.96 

Interface  of  algorithms  at  x  =  5,  x  =  15  mm 

Dl 

20.1875 

0% 

28i 

3s 

31s 

1.7 

HI 

20.1875 

0^ 

288 

6s 

34s 

1.55 

2.2  Absorbing  boundary  condition 

The  absorbing  boundary  condition  was  implemented  for  the  case  of  the  TEiO  mode  in  a  rectangular 
waveguide.  The  offline  identification  procedure  consisting  in  expanding  the  dispersive  part  of  the  impulse 
response  into  series  of  the  Leguerre  polynomials  was  used  do  model  the  dispersive  part  of  the  waveguide. 
The  expansion  can  be  written  in  the  following  form 

p 

hi(t)  =  J2^iLi(t)  (74) 

•  si 

where  Li(t)  denotes  the  Leguerre  polynomial  of  the  order  i.  The  Leguerre  polynomials  have  the  advantage 
of  yielding  an  always  stable  model.  Prior  to  simulation  a  wide  band  signal  shown  in  hg.8  was  propagated 
using  the  PEE  algorithm  for  600  time  steps  and  its  samples  taken  at  two  adjacent  slices  were  recorded. 
Based  on  this  data  the  expansion  coefficients  in  (74)  were  found  and  used  during  the  actual  simulations. 
Fig. 9  shows  the  return  loss  of  the  ABC  in  a  i0.16  by  22.86mm  rectangular  waveguide.  The  solid,  dashed 
and  dashed  dotted  lines  correspond  to  p  =  2,  p  =  5  and  p  =  10  in  expansion  (74),  respectively.  It  is  seen 
that  even  for  the  lowest  approximation  order  the  results  are  excellent  as  the  return  loss  is  lower  than  -80 
dB  over  a  very  wide  frequency  band.  The  worst  performance  is  seen  near  the  cutoff  frequency.  Fig.  10 
shows  the  close  up  of  the  cutoff  region  it  is  seen  that  increasing  the  approximation  order  the  the  quality 
of  the  ABC  improves  quickly.  Since  the  practical  structures  do  not  operate  at  cutoff  frequency  -20  dB 
return  loss  at  this  point  will  have  no  influence  on  actual  modeling.  It  should  be  noted  that  the  numerical 
cost  of  implementing  the  proposed  ABC  is  extremely  low.  The  overhead  is  marginal  as  the  calculation 
of  response  requires  only  a  few  times  the  number  of  expansion  terms  of  extra  floating  point  operations 
and  the  storing  the  a  few  expansion  coefficients.  It  has  to  be  noted  that  this  cost  is  a  few  orders  of 
magnitude  lower  than  in  the  most  recently  published  results  [23,  25, 24,  26]  of  the  other  authors,  where 
form  few  hundreds  to  few  thousand  of  extra  floating  point  operations  at  each  time  step  and  a  similar 
number  of  additional  memory  locations  was  used  to  achieve  a  similar  performance  of  ABC. 
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Figure  8:  The  excitation  used  for  off  line  identification  of  the  ABC 
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Figure  9:  The  return  loes  of  the  ABC  for  a  TEio  mode  in  a  10.16  by  22.86mm  rectangular  waveguide. 
The  solid,  dashed  and  dashed  dotted  lines  correspond  to  two,  five  and  ten  Leguerre  polynomials  used  in 
expansion  (74),  respectively 
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Figure  10:  The  return  loss  of  the  ABC  for  a  TEio  mode  in  a  10.16  by  22.86mm  rectangular  waveguide. 
The  solid,  dashed  and  dashed  dotted  lines  correspond  to  two,  live  and  ten  Leguerre  polynomials  used  in 
expansion  (74),  respectively 
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3  DISCUSSION 

3.1  Algorithm  limitations  and  their  applicability  to  complex  problems 

The  time  domain  algorithms  based  on  the  method  of  moments  and  the  eigenfunctions  expansion  are  in 
general  less  versatile  than  classical  FDTD  schemes.  One  limiting  factor  is  that  the  structure  has  to  be 
entirely  (for  the  TEE)  or  partially  (for  the  PEE)  closed  (  As  discussed  further,  this  restriction  can  be 
overcome  in  the  PEE  if  the  spherical  coordinate  system  is  used).  Also  because  the  eigenfunction  are 
available  in  analytical  form  only  for  simple  domains,  the  TEE  and  PEE  are  restricted  to  the  analysis 
of  the  problems  which  can  be  described  within  such  domains.  Additionally,  the  method  of  moments 
procedure  requires  computation  of  inner  products.  The  numerical  cost  involved  in  this  process  limits  the 
applicability  of  TEE  and  PEE  to  weak  rmd  moderate  inhomogeneities. 

These  drawbacks  can  be  alleviated  in  two  ways,  one  is  to  replace  the  eigenfunction  expansion  with 
the  basis  functions  with  local  support.  This  approach  leads  to  a  formulation  which  may  be  regarded 
time  domain  version  of  the  finite  element  method.  The  other  possibility  is  to  combine  the  PEE  with  a 
classical  FDTD.  The  hybrid  technique  will  be  particularly  useful  in  structures  where  the  inhomogeneities 
are  separated  one  from  the  other  by  homogeneous  regions.  One  problem  with  a  classical  FDTD  method 
is  that  the  field  from  one  inhomogeneity  to  the  other  has  to  be  propagated,  because  of  the  stability 
requirements  the  time  step  size  is  usually  small  and  if  the  smallest  feature  size  is  also  small  then  the 
propagation  of  a  field  through  homogeneous  regions  requires  vast  amounts  of  computer  store  and  time. 
If  the  homogeneous  regions  can  be  treated  with  scalar  PEE  algorithm  great  savings  in  computer  CPU 
time  and  memory  resources  can  be  achieved  as  it  often  suffices  to  propagate  a  few  low  order  modes  and 
additionally  the  scalar  PEE  algorithm  offers  near  ultimate  speed.  The  scalar  PEE  may  be  regarded  here 
as  a  fast  field  propagator.  As  shown  in  the  numerical  tests  the  hybrid  algorithms  are  also  extremely  useful 
to  implement  the  ABC  in  homogeneous  guides  so  their  implementation  can  offer  a  significant  speed  up 
in  the  analysis  of  complex  microwave  structures. 

Examples  of  practical  structures.  It  should  be  noted  that  there  many  practical  complex  structures 
whose  geometries  are  suitable  for  the  application  of  scalar  hybrid  PEE-FDTD.  For  instance  most  compo¬ 
nents  using  rectangular  waveguide  such  as  waveguide  to  coax  transition,  all  kinds  of  E-plane  ad  dielectric 
filters,  septa,  T-  and  E-H  junctions  which  lead  to  3D  problems  can  easily  be  represented  in  a  combination 
of  PEE  and  FDTD  meshes  with  2D  homogeneous  slices  filling  most  of  the  volume.  In  such  structures 
the  speed  and  memory  improvement  should  be  significant.  Also  many  2D  structures,  including  most  of 
planar  transmission  lines  can  be  described  in  a  hybrid  mesh  with  ID  "slices”  located  in  homogeneous 
part  and  FDTD  mesh  located  only  in  the  vicinity  of  strips,  slots  etc.  The  hybrid  approach  solves  here 
the  problem  encountered  in  a  classical  FDTD  where  the  fine  mesh  density  required  in  near  the  small 
geometrical  features  has  also  to  cover  large  homogeneous  regions. 

3.2  FVee  space  problems 

Hybrid  PEE-FDTD  techniques  seem  promising  for  problems  of  field  scattering  or  radiation  in  the  free 
space.  Here  the  best  approach  is  to  place  an  object  in  spherical  domain  and  use  the  FDTD  mesh  inside 
the  sphere  and  the  eigenfunction  expansion  outside.  In  the  PEE  peut  of  the  algorithm  we  discretize  the 
space  along  the  r  direction  rmd  use  spherical  harmonics  for  the  expansion  in  the  ip  and  Q  directions.  In 
spherical  coordinates  an  open  space  may  be  treated  as  a  spherical  waveguide  [14]  which  propagates  the 
TE  and  TM  modes.  For  an  open  space  the  wave  impedance  has  the  resistive  character  only  if  *r  >  n 
where  n  is  the  order  of  the  spherical  Bessel  functions.  The  modes  with  the  resistive  impedance  radiate 
while  those  with  the  imaginary  character  of  the  impedance  are  confined  to  the  region  in  the  neighborhood 
of  the  sources.  Hence  low  order  spherical  harmonics  can  be  used  to  describe  radiation  (scattering)  from 
the  sources  contained  in  a  sphere  of  the  radius  comparable  with  the  shortest  wavelength  of  interest  (ie 
satisfying  fcr  >  n).  Note  that  the  field  expansion  in  the  spherical  coordinates  uses  Fourier  series  in  the 
ip  direction  so  the  FDTD  and  PEE  can  be  easily  interfaced  through  the  FFT.  There  could  be  a  few 
advantages  of  this  approach.  First,  we  can  expect  that  a  hybrid  algorithm  would  require  much  fewer 
unknowns  than  in  the  FDTD.  Second,  ABC’s  can  be  easier  to  implement  on  mode  by  mode  basis  as 
discussed  in  Section  1.7.  A  very  good  performance  can  be  expected.  Third,  the  PEE  can  be  used  as  a 
field  propagator  to  allow  analysis  of  fields  at  larger  distances. 
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3.3  Future  developments 

Having  created  the  theoretical  basis  for  the  algorithms  and  validated  the  key  techniques,  it  is  possible  do 
define  future  research  directions.  There  are  four  topics  which  seems  worth  pursuing. 

1.  A  full  implementation  of  a  scalar  PEE^FDTD  algorithm.  As  shown  this  algorthm  offers  near 
ultimate  speed  requiring  little  computer  store. 

2.  Optimisation  of  the  new  ABC  for  highly  dispersive  structures.  The  results  obtained  so  far  represents 
the  qualitative  and  quantitative  leap  in  speed,  memory  and  performance  compared  to  the  solutions 
currently  used.  The  emphasis  in  future  research  should  be  placed  on  the  reduction  of  the  reflection 
at  cutoff  frequency. 

3.  The  development  of  the  hybrid  PEE-FDTD  scheme  for  open  space.  This  can  expand  the  applica¬ 
bility  of  the  methods  of  moments  in  time  domain. 

4.  A  basic  theoretical  research  in  order  to  develop  a  suitable  set  of  entire  domain  expansion  func¬ 
tions  for  20  inhomogeneous  regions.  The  PEE  algorithm  currently  requires  homogeneous  slices 
in  order  to  achieve  great  speed.  With  a  suitable  basis  constructed  basis,  taking  into  account  the 
inhomogeneities  of  2D  slices  abd  allowing  the  arbitrary  shapes  of  the  region,  the  hybrid  PEE-FDTD 
algorithm  should  give  the  same  versatility  as  FDTD  but  at  a  much  greater  speed  for  practically  all 
microwave  structures  which  can  be  analyzed  by  current  time  domain  methods. 

All  four  topics  are  promissing  from  the  point  of  view  of  numerical  cost  reduction  and  can  be  investigated 
at  the  Technical  University  of  Gdansk. 
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SUMMARY 

Series  function  expansion  algorithms  of  the  tune-domain  analysis  of  boundary  value  problems  are  discussed. 
Electromagnetic  fields  inside  a  structure  under  investigation  are  expanded  into  series  of  basis  functions  and 
the  expansion  coefficients  are  found  by  means  of  the  Galerkin  method.  The  numerical  cost  of  algorithms 
is  discussed  and  a  cost  efficient  approach  is  proposed  for  formulations  using  sine  and  cosine  expansion 
functions.  Compared  with  conventional  time-domain  methods  the  algorithms  described  show  the  time  evol¬ 
ution  of  the  expansion  coefficients  rather  than  the  samples  of  a  physical  continuum  at  discrete  nodes. 


1.  INTRODUCTION 

The  Fourier  transform  introduces  the  equivalence  between  time  and  frequency.  Therefore  each 
problem  involving  one  of  these  quantities  has  two  alternative  formulations.  For  many  years  the 
frequency-domain  formulations  have  been  preferred  for  the  analysis  of  Maxwell  equations.  A 
plethora  of  frequency-domain  numerical  methods  have  been  developed  and  applied  to  various 
boundary  value  problems  in  electromagnetics.  In  recent  years  time-domain  techniques  have 
increasingly  been  gaining  an  audience.  Most  research  in  the  area  of  time-domain  algorithms  for 
the  modelling  of  electromagnetic  fields  was  concerned  with  three  versatile  and  simple  methods. 
These  three  algorithms  are  known  as  the  transmission  line  matrix  (TLM),  spatial  network  (SN) 
and  finite  difference,  time-domain  (FDTD)  methods.  All  three  techniques  are  well  described  in 
the  open  literature  and  therefore  we  cite  here  only  a  few  sources.'*^  In  the  TLM  and  SN  methods 
the  structure  under  investigation  is  treated  as  a  spatial  network  of  transmission  lines  and  the  wave 
propagation  is  described  by  incident  and  reflected  voltage  impulses  on  the  mesh  lines.  In  the 
FDTD  method  the  Maxwell’s  equations  are  discretized  both  in  space  and  time,  the  derivatives 
are  computed  by  means  of  central  difference  scheme  and  the  fields  are  evaluated  at  discrete 
nodes.  Although  FDTD  and  TLM  (SN)  methods  use  different  concepts  and  are  concerned  with 
different  physical  quantities,  a  recent  study  by  Celuch-Marcysiak  and  Gwarek^  shows  that  each 
of  them  can  be  obtained  from  the  other  by  a  sequence  of  suitable  transformations.  Hence,  time- 
domain  methods  currently  used  in  practice  are  formally  equivalent.  Their  salient  feature  is  that 
three-dimensional  space  is  discretized  and  a  mesh  is  formed.  Samples  of  relevant  quantities  at 
mesh  points  are  used  to  represent  a  physical  continuum.  Obviously,  while  solving  complex 
problems  we  have  always  to  contend  with  approximate  answers,  but  nevertheless  sampling  of  the 
solution  in  space  is  only  one,  not  necessarily  the  most  efficient,  form  of  a  discrete  representation 
of  a  physical  continuum.  For  instance  in  the  frequency-domain  methods  this  form  of  approximation 
is  present  in  the  finite  difference  technique*  and  to  a  certain  extent  also  in  the  method  of  lines.* 
Offier  algorithms  use  different  representations  of  fields.  One  of  the  most  popular  frequency- 
domain  techniques  consists  in  field  expansion  in  the  series  of  basis  functions.  Depending  on  the 
algorithm,  expansion  functions  are  defined  locally  or  over  the  entire  domain.  They  can  be  just 
mathematical  functions  or  have  a  physical  mei^ng  of,  for  instance,  modes.  The  expansion 
coefficients  of  are  most  often  found  using  the  methods  known  from  functional  analysis,  such  as 
the  method  of  moments  (Galerkin,  Ritz),  or  least  squares.'*  Expansion  of  the  fields  in  series  of 
basis  fiinctions  underlies  tudi  techniques  as  finite  elements,  coupled  modes,  mode-matching, 
point-matching  and  iterative  eigenfunction  expansion,'**'*  to  mention  only  the  most  poweiful. 
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To  the  author's  best  knowledge,  the  only  attempt  so  far  to  apply  function  expansion  in  the  time- 
domain  analysis  of  shielded  microwave  structures  was  published  by  Nam  et  al7  The  expansion 
function  used  in  that  case  was  obtained  using  frequency-dependent  equations  and  the  algorithm 
was  far  more  complicated  than  that  of  classical  time-domain  methods. 

This  introduction  shows  that,  as  far  as  the  approximation  techniques  are  concerned,  the  time- 
domain  methods  of  the  modelling  of  electromagnetic  fields  have  used  so  far  only  a  small  fragment 
of  a  range  of  choices  existing  in  the  frequency  domain.  In  view  of  that  fact  the  aim  of  this  paper 
is  to  look  at  time-domain  algorithms  from  the  point  of  view  of  the  alternative  representation  of 
fields,  in  order  to  find  formulations  which  may  broaden  the  range  of  options  available  for  the 
time-domain  analysis  of  complex  problems  of  electromagnetics. 


2.  ANALYSIS 


We  shall  be  concerned  with  the  time-domain  analysis  of  shielded  structures.  Of  particular  interest 
are  resonators  and  cylindrical  waveguides  inhomogeneously  loaded  with  non-dispersive  isotropic 
materials.  One  possible  structure  is  shown  in  Hgure  1.  It  is  assumed  that  the  bounding  walls  are 
perfect  electric  or  magnetic  conductors  and  that  the  inhomogeneity  is  described  by  the  relative 
primitivity,  c„  .md  permeability,  p,,  both  being  in  general  functions  of  all  three  space  coordinates. 
Under  these  assumptions  the  Maxwell’s  equations  are  given  by 


i-E  =  — 5— 
at  tffir{x,y,z) 


V  X  H 


- ? - 


where  co  and  are  the  permittivity  and  permeability  of  the  free  space. 
These  equations  can  be  written  using  the  following  abbreviated  notation 


(1) 


(2) 


where  if,  =  V  x  (•),  =  ~^fv^ix,y,z)  V  x  and  /  and  g  are  vector 

functions  representing  the  electric  and  magnetic  field,  respectively.  The  above  notation  will  be 
used  henceforth. 

Replacing  the  time  derivatives  with  a  finite  difference  equivalents  we  get  the  following  time¬ 
matching  equations 


/"=/*•-> -H  Arif, g-- >'2 

=  Arifi/"  (3) 

The  unknown  functions  /,  g  are  now  expanded  into  series  of  functions. 


Figure  1.  A  powible  geooieuy  of  the  probkin 
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/"  ”  2  /XJC.y.^)  «"  *  2  l>7g^x,y,z)  (4) 

The  set  of  expansion  functions  is  assumed  to  be  complete  and  the  functions  are  linearly  indepen¬ 
dent.  The  functions  are  defined  on  the  entire  domain  or  have  local  supports  but  they  are  in 
general  the  time-independent  functions  of  all  three  spatial  variables.  Substituting  (4)  into  (3)  we 
get 

2 -  2  «r7/ + Ar  2 

+  (5) 

Taking  the  inner  product  of  (3)  with  the  expansion  functions  //  and  gt  results  in 

+  Ar('e2/".ft>  (6) 

Using  (S)  the  above  equations  can  be  cast  into  the  following  matrix  form 

where  a  and  b  are  the  vectors  containing  expansion  coefficients  and  A,  B,  C,  D  are  matrices  with 
elements  given  by  the  following  inner  products 

~  ifftfi)  (8) 

All  matrices  defined  above  are  in  general  dense.  If  the  expansion  functions  are  orthonormal  then 
C  and  D  are  identity  matrices  and  (7)  becomes 

a"  =  a"-* -»•  A/ Ab"->'2 

So  far  we  have  presented  the  algorithms  in  which  the  time  is  discretized  and  the  function 
expansion  is  done  in  three  dimensions.  We  shall  call  them  total  expansion  algorithms.  Another 
version  of  the  expansion  algorithm  is  obtained  if  the  discretization  is  in  time  and  one  selected 
coordinate  and  the  expansion  is  done  with  respect  to  two  remaining  coordinates.  This  approach 
combines  the  FDTD  method  with  the  expansion  algorithms  decribed  above.  The  space  b  sliced 
into  subdomains  (Figure  2)  and  the  fields  are  expanded  on  each  subdomain  (slice)  into  series  of 


Rgun  2.  Diviskm  of  the  seonetiy  for  the  pettiai  expamioa  alioriibint 
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expansion  functions.  This  type  of  expansion  we  shall  call  partial.  To  obtain  the  time-marching 
equations  for  such  a  case  we  have  to  introduce  the  finite  difference  scheme  for  the  calculation 
of  the  space  derivatives  with  respect  to  the  discretized  variable.  Suppose  the  structure  was  divided 
into  K  dices  in  the  z-direction  and  the  slices  are  uniformly  spaced  by  the  distance  Ad.  To  enable 
the  application  of  the  central  finite  difference  scheme  for  the  calculation  of  the  derivatives  in  the 
z-direction  we  use  the  technique  used  in  the  FDTD  method.  In  this  technique^  different  field 
components  are  defined  on  two  meshes,  each  off  by  half  the  space  step  from  the  other.  This 
arrangement  is  known  as  Yee’s  mesh.  In  the  present  context  Yee’s  mesh  is  given  in  one  spatial 
dimension,  which  means  that  relevant  field  components  are  defined  on  the  slices  that  are  off  by 
half  the  space  step.  More  precisely  the  £^,  £^, are  given  at  z  -  kAd  and  //f, //*,£*  are 
defined  for  z  =  (k  +  l/2)Ad.  Bearing  this  convention  in  mind  we  expand  fields  at  a  suitably 
situated  slice  according  to 

/*  =  2  flr*4(x,y)  «*  =  2  i . ( lo) 

The  derivatives  in  the  z-direction  are  approximated  by 


Ad 


(11) 


Calculation  of  the  z-derivatives  involves  operations  on  the  functions  defined  on  the  adjacent 
slices.  Therefore,  instead  of  equations  (2)  we  get 


£**  =  nA  +  iei(/* -/*-,)  (12) 


Where  we  have  split  operators  and  X2  into  transverse  and  z-  part  according  to 


J  \V,X(-) 
c<,e*(z,y) 

( 

- 

«o€*(z,y)Ad  * 
-1 

mi4*(x,y)Ad 


i,  X  (•) 


In  the  above  equations  by  i«  we  have  denoted  a  unit  vector  in  the  z-direction. 

Replacing  time  derivatives  by  the  finite  difference  scheme,  expanding  functions  in  (12)  and 
taking  the  inner  products  for  each  slice  we  arrive  at  the  set  of  equations  similar  to  (7)  with  the 
block-diagonal  coefficient  matrices  given  by 


A  =  qdiag  [A'*,A"*)  B  =  qdiag 

C  =  qdiag  (C*J  D  =  qdiag  (D*]  (14) 


The  elements  of  the  !>ubmatrices  are  given 


Ai* » (cn  -  Ai*  - 

Bv*  -  ((%  +  ^)4,g,.)  By  -  - 

CJ”  ()>*./<*) 

(15) 


For  orthonormal  basis  functions  matrices  C  and  O  become  identity  matrices  and  we  get  again 
(9).  Thus  both  total  and  partial  time-domain  expansion  algorithms  are  formally  described  by  the 
same  set  of  equations  (7)  or  (9).  However,  one  important  difference  between  the  total  and  partial 
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expansion  algorithms  is  that  in  the  former  the  matrices  (8)  are  in  general  dense  while  in  the 
latter  the  matrices  (15)  are  always  quasi>diagonal. 

Compared  with  the  FDTD  method  the  expansion  algorithms  show  the  time  evolution  of  the 
expansion  coefficients  rather  then  field  components  at  nodes.  Obviously  the  expansion  algorithms 
have  to  use  timestep  small  enough  to  prevent  the  instability.  At  this  point  it  is  not  possible  to 
give  one  stability  criterion  for  the  algorithms  discussed  as  it  depends  on  the  choice  of  the  expansion 
functions. 


3.  BASIS  FUNCTIONS 

Before  we  discuss  the  numerical  cost  involved  in  the  time  domain  expansion  algorithms  we  shall 
consider  the  choice  of  expansion  functions.  As  we  indicated  earlier,  in  both  total  and  partial 
expansion  algorithms  individual  basis  fiinctions  can  have  local  support  or  be  defined  in  the  entire 
domain  or  subdomain  (slice).  The  advantage  of  the  first  approach  is  that  it  increases  the  qnrsity 
of  the  matrices.  On  the  other  hand  choosing  the  basis  functions  which  have  the  space  distribution 
of  the  relevant  field  components  of  modes  in  the  corresponding  homogeneous  problem,  we  obtain 
the  following  advantages 

(1)  expansion  functions  satisfy  the  boundary  conditions 

(2)  owing  to  the  orthogonality  of  mode  functions  in  homogeneous  problems,  the  inner  products 
at  homogeneous  slices  (Figure  2)  in  the  partial  expansion  algorithms  result  in  diagonal  sub¬ 
matrices  A'*,  A''*.  B"*,  B'** 

(3)  the  time-marching  algorithms  receive  a  new  physical  interpretation.  They  can  now  be  viewed 
as  the  equations  describing  mode  coupling  due  to  inhomogeneity.  This  interpretation  allows 
a  straightforward  investigation  of  mode  interaction  in  loaded  guides  and  resonators 

(4)  for  small  inhomogeneities  only  a  few  expansion  functions  will  be  sufficient  to  approximate 
field  distribution  with  good  accuracy 

A  suitable  choice  for  the  time-independent  basis  functions  which  have  the  distribution  of  modes 
of  a  homogeneous  structure  are  the  eigenfiinctions  of  the  Laplace  operator.  The  detailed  discussion 
of  the  choice  of  the  complete  basis  for  arbitrarily  shaped  domains  and  the  relation  between  the 
eigenfunctions  of  Laplace  operator  and  modes  of  a  homogeneous  structure  is  given  in  Reference 
IS.  This  discussion  shows  that  for  the  total  expansion  algorithms  and  arbitrarily  shaped  region 
closed  by  the  perfect  electric  conductor  5  we  may  use  the  following  representation  of  fields 

E  =  2 

H-2(«»A  +  h;h;)  (16) 

where  Ci  and  h<  are  solenoidal  eigenfunctions  of  a  vector  Laplace  operator 

V^  +  Xe  =  0  V*h  +  Xh  =  0  (17) 

with  suitable  boundary  conditions,  and  ci  and  h!  are  potential  functions  obtained  from  eigen¬ 
functions  of  a  scalar  l^place  operator 


V24»  -Hwh  =  0  (18) 

with  Dirichlet  (for  the  derivation  of  e!)  and  Neumann  (for  h))  boundary  conditions.  Here  functions 
tt  and  li«  correspond  directly  to  the  modes  of  a  homogeneous  resonator.  In  the  case  of  the 
partial  eigenfunction  expansion  a  suitable  form  of  representation  involving  tune-independent  basis 
functions  is 


E  «  2  +  a^e^  + 

H -2  +  + 


(19) 
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Subscripts  i  and  /  refer  to  z  and  transverse  component,  respectively.  Functions  t,„  ba.  t'u  are 
obtained  from  the  eigenfunctions  of  a  scalar  Li^laiix  curator  with  Oirichlet  boundary  conditions 
and  functions  b^i.  b'^  are  derived  from  the  eigenfunctions  of  a  scalar  Laplace  operator  with 
Neumann  bounduy  conations.  In  terms  of  TE  and  TM  modes,  functions  e,i,  b,<,  t'n  contribute 
to  TM  modes  in  a  homogeneous  guide  and  b^i.  Cm,  b^i  produce  TE  fields. 

For  simple  regions,  5,  it  may  be  also  possible  to  expand  each  field  component  separately 
using  the  eigenfunctions  of  a  sc^ar  Laplace  operator  with  mixed  Dirichlet/Neumann  boundary 
conditions. 


4.  NUMERICAL  COSTS  OF  THE  TIME-DOMAIN  EXPANSION  ALGORITHMS 

One  drawback  of  the  expansion  algorithms  presented  in  the  previous  section  is  that  they  may 
lead  to  higher  numerical  cost  than  FDTD  a^  TLM.  The  inner  (voducts  appearing  in  (8)  and 
(IS)  are  independent  of  time  and  can  be  stored  in  look-up  tables.  If  the  cost  of  the  calralation 
of  inner  (»oducts  is  neglected  in  the  algorithms  then  the  cost  of  one  timestep  of  the  algorithms 
is  determined  by  the  cost  of  matrix  multiplicatitm.  Depending  on  the  basis  functions  used  (local 
support  versus  entire  domain)  the  overall  cost  may  vary  considerably.  Generally  speaking,  assuming 
that  expansion  is  done  using  L  eigenfunctions,  the  cost  of  one  timestep  is  of  order  O(L^).  If  the 
matrices  are  sparse  the  cost  is  lower.  In  the  FDTD  and  TLM  method  with  N  nodes,  the  numerical 
cost  is  of  order  0(Af).  Consequently,  expansion  techniques  are  comparable  in  terms  of  numerical 
cost  and  memory  to  classical  time-domain  algorithms  when  ~  N.  This  condition  will  easily  be 
fulfilled  in  slightly  and  moderately  perturbed  homogeneous  structures  when  the  basis  functions 
have  the  space  dstribution  of  the  relevant  field  components  of  modes  in  the  corre^nding 
homogeneous  problem.  A  similar  conclusion  can  be  derived  regarding  memory  requirements. 

At  this  point  it  is  important  to  note  that  a  very  efficient  implementation  of  the  expansion 
algorithms  may  be  obtained  if  the  expansion  functions  are  sine  and  cosines.  Let  us  consider  the 
total  expansion  algorithm  described  by  (9).  Equations  (6)  imply  that  at  each  step  one  evaluates 
the  inner  products  ft)  and  (ifj  /'*•{•)  ainri  these  inner  products  are  used  to  update 

expansion  coefficients.  For  sine  and  cosine  functions  the  inner  product  for  all  testing  functions 
can  be  computed  in  a  very  efficient  way  using  the  technique  described  in  Reference  14.  In  this 
technique  tte  off  inner  pr^ucts  are  computed  in  one  step  in  a  sequence  of  inverse  and  forward 
FFTs.  In  a  nutshell  the  procedure  is 

•  using  invent  3-D  FFT  and  and  a"  calculate  f  and  their  spatial  derivatives 

•  Compute  functions  ifig""*'*  and  ^2 1 

•  using /onvard  3-D  FFT  compute  all  inner  products  (^ig""*'*, /*>  and  (SEjf  ,gi) 

The  numerical  cost  of  such  computations  is  relatively  low  and  there  is  no  need  to  create  look¬ 
up  matrices  A  and  B.  A  result  the  time-domain  algorithms  can  be  implemented  using  very  little 
computer  storage  and  with  a  speed  comparable  to  that  of  the  corresponding  FDTD  or  TLM 
methods. 

It  is  beyond  the  scope  of  this  paper  to  discuss  in  detail  the  numerical  implementations  of  the 
algorithms  presented  herein.  However,  both  total  and  partial  expansion  algorithms  have  been 
implemented  and  tested  by  the  author  for  case  of  the  rectangular  resonators  containing  dielectric 
slabs.  The  eigenfunaions  of  the  Laplace  operator  with  Dirichlet  or  Neumann  boundary  conditions 
were  used  as  the  expansion  functions  and  the  inner  products  were  calculated  using  fast  Fourier 
transform  a^troach  described  in  the  above  section.  The  algorithms  have  been  found  stable.  Some 
numerical  re^ts  for  the  expansion  algorithms  with  sine  and  cosine  expansion  functions  are  given 
in  Table  1  and  Reference  16.  Table  1  compares  the  results  obtained  for  fundamental  mode  in  a 
rectangular  resonator  shown  in  Figure  3  using  the  total  and  partial  expansions  algorithms  with 
the  data  FDTD  published  in  Reference  S  and  rigorous  methods.  It  can  be  noted  that  the  results 
in  the  ubie  obl^ned  with  the  expanskm  algorithms  are  in  a  very  good  agreement  with  the 
reference  values  which  confirm  the  validity  of  the  analysis  proposed  in  this  paper. 


5.  CONCLUSIONS 

Algorithms  of  the  time-domain  analysis  of  inhomogeneously  loaded  microwave  structures  have 
been  described.  The  methods  propoi^  are  based  on  the  expansion  of  fields  into  complete  series 
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Tabic  I.  Comparison  of  the  results  for  the  normalized  resonant  frequencies  A//X  obtained 
using  expansion  algorithms  with  the  published  data’  for  the  rectangular  resonator  shown 

in  Figure  3 


TRM 
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expansion 
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expansion 

FDTD’ 

w  »  0 

0-07S11  (exact) 

0-07502 

0-07511 

0-0750 

w  -  2A/ 

0-0322 

0-05253 

0-05153 

0-0517 

w  w  4A/ 

0-0445 

0-04456 

0-04433 

0-0442 

Figure  3.  The  geometry  of  the  problem  used  in  the  numerical  test  (a  -  12AI.  b  •  6^.  c  >  8d/,  c,  3-75) 


of  basis  functions.  The  resulting  equations  show  the  time  evolution  of  the  expansion  coefficients. 
With  a  suitable  choice  of  basis  functions  this  feature  allows  one  to  investigate  propagation  of 
separate  modes  and  their  mutual  interactions. 
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ABSTRACT 

The  paper  presents  the  derivation  of  the  stability 
condition  for  the  total  rigenfhnction  expansion  al¬ 
gorithm  for  the  solution  of  Maxwell’s  equations. 
The  algorithm  uses  the  method  of  moment  solu¬ 
tion  for  space  and  leap-frog  differencing  scheme  for 
time.  The  stability  condition  of  the  algorithm  is 
derived  by  investigating  the  properties  of  operators 
in  suitably  defined  Hilbert  spaces.  The  approach  is 
general  and  can  be  used  for  other  iterative  algori¬ 
thms. 

1.  INTRODUCTION 

At  present  there  are  two  algorithms,  namely  the 
TLM  and  FD-TD  [1],[2],  wluch  dominate  the  arena 
of  time  domain  solution  of  Maxwell’s  equations.  In 
both  methods  the  space  is  discretized  and  derivati¬ 
ves  are  approximated  by  finite  differences.  H,  howe¬ 
ver,  the  sampling  of  at  the  me'h  nodes  is  replaced 
by  the  series  expansion  of  fields  and  the  expansion 
confidents  are  found  using  the  method  of  moment, 
new  algorithms  can  be  created.  Two  such  algori¬ 
thms  have  been  recently  proposed  by  the  author  of 
this  contribution  [3],  [4]. 

A  simplest  time  domain  algorithm  based  on  the 
concept  briefly  described  above  is  called  a  total 
dgenfunction  expansion  (TEE)  [3],  [4].  Let  ns  con¬ 
sider  a  set  of  coupled  differential  equations  reflec¬ 
ting  the  form  of  Maxwdl  equatkms 


-  Ug 

§^9  =  Uf  (1) 

^  ^  ^  {•) 
and  /  and  g  are  vector  functions  representing  the 
electric  and  magnetic  field,  respectively.  Instead  of 
discretizing  the  fields,  equations  (1)  are  solved  by 
means  of  the  method  of  moments  and  only  time 
derivatives  are  calculated  using  the  finite  differen¬ 
cing  scheme.  Repladng  the  time  derivatives  with 
finite  difference  equivalents  we  get  the  following 
time  marching  equations 

/"  =  /"-* -H 

gn+i/i  ^  +  (2) 

The  unknown  functions  f,g  are  now  expanded  into 
series  of  functions. 

(3) 

Expansion  function  are  defined  on  the  entire 
domain.  A  sensible  choice  for  the  electromagne¬ 
tic  fields  are  the  rigenfunctions  of  Laplace  opera¬ 
tor  with  suitable  boundary  conditions.  Applying  a 
standard  method  of  moment  procedure  by  taking 
the  inner  product  of  (2)  with  the  expansion  func¬ 
tions  fi  and  gi  we  get 
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<pji>  =  <r~'ji> 

+  Al<L3r,Si>(4) 

The  above  equations  can  be  cast  into  the  fdlowing 
matrix  form 

ft"  =  ft"-*  + 

jn+i/3  ^  fc»-»/a  + 

where  ft  and  ft  are  the  vectors  containing  expansion 
coefficients  and  ^  •  £  <  £  i  £  matrices  with  ele¬ 
ments  given  by  the  following  inner  products 

Aij  —  <  ^  Bij  ^ 

=  <  /ji/i  >  ^i]  =<  ftjtfti  ^  (6) 

2.  STABILITY  CONDITION 

The  algorithm  described  above  was  proven  to  pro¬ 
vide  an  accurate  solution  for  Inhomogeneously  lo¬ 
aded  rectangular  resonator  [3]i  [4].  However,  during 
the  numerical  tests  it  was  observed  that  the  maxi¬ 
mal  time  step  size  depends  on  the  number  of  func¬ 
tions  used  in  the  expansion.  This  time  step  size  is 
a  key  factor  which  allows  one  to  derive  the  stabi¬ 
lity  criterion  for  the  algorithm.  For  the  FD-TD  and 
TLM  method  the  stability  criterion  is  know  as  a 
Courant  conditions,  but  the  stability  of  TEE  algo¬ 
rithm  has  never  been  examined  yet.  This  contribu¬ 
tion  will  present  the  results  of  the  stability  analysis 
of  the  TEE  algorithm 

Presently  we  shall  derive  the  stability  criterion 
for  the  the  total  expansion  using  trygonometric  se¬ 
ries  expansion.  Before  deriving  the  stability  crite¬ 
rion  let  us  present  problem  (1)  in  the  form  of  a 
hyperbolic  equation.  Thking  the  time  derivative  of 
the  first  equation  and  applying  operator  Li  to  the 
second  equations  we  obtain. 

^/  =  l|Io/  (7) 


Putting  L  =  -LiLj  we  get 


I5/+I./.0 


(8) 


L  is  an  elliptic  differential  operator  with  positive 
coefficients  so  (8)  is  hyperbolic  differential  equation. 
Stability  of  the  differencing  schemes  can  be  studied 
using  the  methods  of  functional  analysis.  In  [5]  such 
an  approach  was  used  for  the  case  of  classical  finite 
difference  representation  of  L.  The  theorems  used 
in  that  case  are  also  applicable  to  this  method.  To 
investigate  the  stability  time  marching  algorithms 
for  the  hyperbolic  equations  it  is  useful  to  present 
a  problem  in  a  canonical  form 


(I-|-A’fR)|j/-HA/  =  0 


(9) 


Where  I  is  the  identity  operator. 

The  time  marching  algorithm  for  the  above  pro¬ 
blem  is  stable  if  the  foUowing  conditions  are  fulfil¬ 
led  [5]: 

A  =  A’  >  0,  R  =  R*  >  0  (10) 

R  -  0.5A  >  0  (11) 

In  other  words  for  the  time  marching  adgorithm  to 
be  stable  it  is  sufficient  that  both  operators  A  and 
R  are  self  adjoint  and  positive  and  additionally  the 
operator  R  —  0.5A  is  nonnegative.  The  canonical  | 
form  (9)  is  obtained  from  (8)  by  simply  multiplying  i 
it  by  2  and  writing  the  result  as 

Comparing  (12)  with  (9)  we  get  R  =  I/Ah  and 
A  =  2L  It  can  readily  be  verified  that  operator  L 
is  symmetric  and  positive.  It  suffices  to  verify  the 
condition  (11).  This  condition  is  fulfilled  when 


fell  ^ 


or 


At  < 


IILII 


(13)1 


(H) 
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I  Men  that  the  stability  depends  on  the  norm  of 
t  operator  L.  The  norm  of  the  operator  depends 
I  the  space  it  acts  in.  Let  ns  estimate  the  norm  of 
I  operator  L  in  a  finite  dimensional  space 
ned  by  trigonometric  functions. 

Suppose  that  the  inhomogenoty  is  located  in- 
de  a  cubical  resonator  with  the  dimondnn#  txtxi. 
ke  expansion  functions  are  normalized  products 
f  sines  and  cosines  of  the  form 

,  ixf  kit(  ....  .  . 

sm  -j-  or  cos  -pi  t,ifc  <  Nu  (15) 

Let  P  denote  the  domain  of  the  ori^al  problem 
given  by  equations  (1).  The  basis  functions  (15) 
ipan  a  finite  dimensional  space  Hnh  C  P  in  which 
the  approximate  solution  is  sought  for.  In  such  con¬ 
structed  space  we  calculate  the  upper  bound  of  the 
norm  of  operator  ||L||.  Note  that  ||L||  <  ||Lnill- 
Where 


)  V  X  V  X  (•)  (16) 


[ftttd 


<m<n  =  inf<,(x,y,z), 

Mmi„  =  inf/ir(z,y,a:)  x.y.sgO  (17) 

The  norm  of  operator  L  in  Hftff  can  now  readily 
be  found 


Ill'll  <l|I.m||<r, 


mo* 


(18) 


where  «««  =  (€oPo«minMiiMn)"*^’  w  th«  maximal 
vdodty  which  a  plane  wave  attains  in  the  medium 
filling  the  structure.  Using  the  above  estimation  we 
get  the  following  stability  condition 

At  < - - - p  (19) 

U  a  n  is  a  a  perpendicular  paraUelpied  with  the  di¬ 
mensions  a  x  6  x  /  smd  the  upper  bound  for  t',/,min 
the  trigonometric  expansion  functions  is  Kmi 
Njii  then  the  condition  (19)  becomes 

At  < - \  ..  (20) 


3.  CONCLUSIONS 

The  derivation  of  the  stability  condition  for  the 
total  eigenfunction  expansion  algorithm  for  the  so¬ 
lution  of  Maxwell’s  equations  is  presented..  The  sta¬ 
bility  condition  of  the  algorithm  is  derived  using 
the  methods  of  functional  analysis.  The  approach 
is  general  and  can  be  used  for  other  iterative  algo¬ 
rithms. 
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ABSTRACT 

This  contribution  presents  the  derivation  of  the  stability  condition  for  time 
domain  method  of  moment  for  linear  hyperbolic  differential  equations  The 
algorithm  uses  the  method  of  moment  approach  for  space  variables  and  finite 
difference  scheme  for  time.  The  stability  condition  of  the  algorithm  is  derived 
by  investigating  the  properties  of  operators  in  suitably  defined  Hilbert  spaces. 
The  irethod  discussed  in  this  contribution  has  been  developed  in  context  of 
the  time  domain  analysis  of  Maxwell’s  equations  but  the  approach  is  general 
and  can  be  used  for  other  explicit  algorithms. 

INTRODUCTION 

Functional  analysis  is  one  of  the  most  powerful  tools  of  the  theoretical 
investigation  of  the  basic  properties  of  numerical  methods.  The  methods  of 
functional  analysis  are  commonly  used  in  the  mathematical  physics,  numer¬ 
ical  mathematics  and  computer  science  but  seldom  in  engineering.  At  the 
same  time  the  engineering  creates  demand  for  new  more  efficient  numerical 
methods  which  would  provide  a  sufficiently  accurate  solution  as  fast  as  pos¬ 
sible.  This  results  in  the  constant  improvements  of  published  algorithms  by 
researchers  who  adapt  them  to  their  particular  needs  without  investigation 
of  the  properties  of  modified  algorithms.  For  instance,  explicit  algorithms 
for  the  solution  of  initial  value  problems  have  recently  received  mudi  atten¬ 
tion  among  researchers  involved  in  the  numerical  analysis  of  electromagnetic 
fields.  Two  methods  belonging  to  this  class,  known  as  finite  difference-time 


clonuiin  (FDTD)  and  transmission  line  matrix  (TLM)  algorithms  have  in¬ 
tensively  been  developed  in  the  last  decade.  Their  salient  feature  is  that 
electromagnetic  field  is  analyzed  in  the  time  domain  and  the  samples  of  rel¬ 
evant  physical  quantities  at  nodes  loavted  at  the  discrete  points  in  space  are 
used  to  represent  a  physictd  continuum.  Tliese  two  methods  are  constantly 
being  improved.  The  improvements  include  the  application  of  graded  meshes 
or  non  orthogonal  cells,  application  of  local  approximations  or  extension  of 
the  basic  algorithms  to  the  new  class  of  materitUs  such  as  ferrites  or  disper¬ 
sive  media.  Also  new  concepts  of  space  representation  of  fields  have  been 
introduced.  The  sampling  at  discrete  points  can  be  replaced  by  the  expan¬ 
sion  into  the  series  of  basis  functions  and  the  expansion  coefficients  found  by 
the  method  of  moments  procedure. 

Recognizing  the  progress  achieved  in  the  recent  years  in  the  time  domain 
analysis  of  electromagnetic  fields,  it  should  be  noted,  that  the  explicit  al¬ 
gorithms  underlaying  these  methods  are  not  unconditionally  stable  and  the 
improvements  introduced  to  algoritluns  affect  their  stability.  In  tliis  contri¬ 
bution  we  shall  present  how  the  eflfects  of  the  algorithm  modifications  can  be 
investigated  using  the  functional  analysis. 

STABILITY  ANALYSIS  OF  EXPLICIT  TIME  DOMAIN 

ALGORITHMS 


Let  us  consider  a  hyperbolic  differential  equation 


(1) 


where  L  is  a  linear  elliptic  diflferential  operator  with  positive  coefficients.  The 
hyperbolic  equation  of  this  type,  supplemented  by  conditions  at  t  =  0  can 
be  solved  for  t  >  0  using  a  classical  finite  difference  explicit  algorithm.  For 
L  being  a  Laplacian  the  stability  criterion  for  the  algorithm  is  known  as  the 
Courant  condition.  For  other  operators  it  is  convenient  to  use  the  methods 
of  functional  analysis.  In  [2|  sudi  an  approach  was  used  for  the  case  of  finite 
difference  representation  of  L.  The  theorems  used  in  that  case  are  general 
so  it  is  very  instructive  to  show  how  they  can  be  applied  for  other  expUcit 
algorithms. 


To  investigate  the  stability  of  a  time  marching  algorithms  for  the  hyper¬ 
bolic  equations  it  is  useful  to  present  a  problem  in  a  canonical  form 


(I  +  AnR)^/  +  A/  =  0 


(2) 


Where  I  is  the  identity  operator. 

The  time  marching  algorithm  for  the  above  problem  is  stable  if  the  following 
conditions  are  fulfilled  [2): 


A  =  A*  >  0,  R  =  R*  >  0 
R~0.5A>0 


(3) 

(4) 


In  other  words  for  the  explicit  algorithm  to  be  stable  it  is  sufficient  that  both 
operators  A  and  R  be  self  adjoint  and  positive  and  additionally  the  operator 
R  —  0.5A  be  nonnegative. 

A  linear  operator  F  defined  in  a  Hilbert  space  (//,  <  •  >)  is  self  adjoint  if 

for  any  x,y  &  H 

<  Fx,  y  >=<  X,  Fy  >*  (5) 

An  operator  F  is  positive  P  >  0  (or  nonnegative  F  >  0  )  when  for  all 
X  €  //,  X  #  0  we  have 


<  Fx,x  >  >  0  or  <  Fx,  x  >  >  0 


(6) 


The  canonical  form  (2)  is  obtained  from  (1)  by  simply  multiplying  it  by  2 
and  writing  the  result  as 


6.H 

Comparing  (7)  with  (2)  we  get  R  =  I/A*t  and  A  =  2L 

If  operator  L  is  symmetric  and  positive  than  the  stability  condition  is 

ll^ll  >  IWI 


or 


At  < 


Tin 


(7) 


(8) 


(9) 


It  is  seen  that  the  stability  depends  on  the  norm  of  the  operator  L.  The 
norm  of  the  operator  depends  on  the  space  it  acts  in. 


STABILITY  ANALYSIS  FOR  THE  TIME  DOMAIN  METHOD 

OF  MOMENTS 


Let  us  consider  a  one  dimensional  second  order  equation 

=  (10) 

fix,  to)  =  foix),  fix  =  0)  =  /(z  =  a)  =  0  (11) 

where  6(z)  >  0  is  a  time  independent  continuous  function  of  z,  Instead  of 
using  the  finite  difference  representation  of  ^  let  us  combine  the  explicit 
algorithm  with  the  method  of  moments.  To  this  end  we  will  use  the  finite 
differences  for  the  approximation  of  time  derivatives,  expand  the  fimction 
fix)  into  series  of  sines 


/(^)  =  ^CiSiniixx/a) 


and  use  the  inner  product 


u,  V  >=  /  uvdx 
Jo 


to  find  the  expansion  coefficient  at  any  instance  of  time.  (A  detailed  deriva¬ 
tion  of  the  time  domain  method  of  moments  for  Maxwell’s  equations  can  be 
foimd  in  (IJ) 

It  can  easily  be  verified  that  operator 

L  =  (14) 

is  positive  and  self  adjoint.  This  case  was  considered  previously  so  we  may 
conclude  that  the  algoritlim  is  stable  if 


At  < 


7i^ 


At  this  point  it  is  necessary  to  estimate  the  norm  of  L.  The  problem  is 
defined  in  the  Hilbert  space  spanned  over  sine  functions.  The  norm  of  L  in 
such  a  space  can  be  estimated  as  follows 


l|L||  <  ULmII  = 


<  Lmx.x  > 


^max  o 

or 


(16) 


where  bmai  is  the  maximal  value  of  6(x)  over  the  interval  0  <  x  <  a. 

We  may  conclude  that  the  explicit  algorithm  combined  with  the  method  of 
moment  with  sine  series  will  be  stable  if  the  time  step  is  chosen  such  that 


At  < 


ITT 


Vbmax 


(17) 


Note  that  maximal  time  step  is  inversely  proportional  to  number  of  basis 
functions. 

Obviously,  the  same  procedure  can  be  applied  to  other  types  of  expansion 
functions,  including  for  instance  finite  elements.  It  is  important  to  note 
however  that  the  time  step  in  the  explicit  algorithm  depends  not  only  on  the 
operator  (equation)  solved  but  also  on  the  way  the  approximation  of  space 
is  constructed. 


CONCLUSIONS 

The  application  of  the  functional  analysis  to  the  investigation  of  the  stability 
of  time  domain  algorithms  has  been  presented.  The  method  can  easily  be 
applied  to  the  investigation  of  the  properties  of  novel  time  domain  schemes 
for  Maxwell’s  equations  such  as  the  ones  proposed  in  [1). 
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ABSTRACT:  A  new  algorithm  for  the  time  domain  analysis  of  electromagnetic  waves  in  shielded 
structures  is  presented.  The  algorithm  combines  the  FDTD  with  a  recently  developed  partial 
eigenfunction  expansion  (PEE)  scheme  to  obtain  acceleration  in  numerical  calculation  and  savings 
in  computer  memory.  An  example  of  the  application  of  the  algorithm  is  presented  showing  an 
overall  speed  improvement. 


INTRODUCTION 

Time  domain  techniques  are  gaining  increased  popularity  in  the  analysis  of  electromagnetic 
waves  because  of  their  ability  to  treat  complicated  geometries  over  a  wide  frequency  range.  Be¬ 
cause  the  most  popular  time  domain  algorithms  are  explicit,  the  stability  of  the  algorithm  puts  a 
restriction  on  the  maximum  allowable  time  step.  Additionally,  to  obtain  fine  resolution  of  fields 
near  singularities,  mesh  size  is  reduced  thereby  increasing  the  computer  storage.  The  computation 
time  and  memory  requirements  are  therefore  critic2d  parameters  for  time  domain  algorithms.  Be¬ 
cause  of  that  some  research  effort  has  recently  been  devoted  to  the  acceleration  of  the  traditional 
methods  by  the  application  of  the  signal  processing  [2]  techniques,  graded  mesh  schemes  [1]  or 
elimination  of  redundtmt  field  components  [3], 

This  contribution  presents  a  new  algorithm  for  shielded  structures  which  consists  in  the  replace¬ 
ment  of  the  FDTD  calculation  in  homogeneous  shielded  subregions  by  eigenfunction  expansion. 
The  eigenfunction  expansion  schemes  in  time  domain  proposed  in  [4]  rely  on  the  expansion  of 
the  unknown  functions  of  selected  space  variables  into  series  of  betsis  functions  and  application 
of  method  of  moment  procedure  to  find  the  expansion  coefficients.  One  version  of  the  algorithm 
csJled  the  Partial  Eigenfunction  Expansion  (PEE),  is  obtained  if  the  function  expansion  is  done 
with  respect  to  two  selected  space  coordinates  while  the  third  spatial  coordinate  and  time  are 
discretized  in  a  way  analogical  to  the  finite  difference  scheme.  By  combining  such  am  approach 
with  a  classical  FDTD  algorithm  an  improvement  in  speed  can  be  obtained  for  an  important  class 
of  shielded  structures. 


FORMULATION 

As  the  FDTD  is  a  well  known  technique,  let  us  begin  with  the  presentation  of  the  second  scheme, 
namely  the  Partial  Eigenfunction  Expansion  (PEE).  Consider  a  dielectric  inhomogeneity  located  in 
a  rectangular  waveguide  (fig.l).  In  the  PEE  the  computational  space  is  sliced  into  subdomains  and 
the  fields  ue  expanded  on  each  subdomain  (slice)  into  series  of  expansion  functions  which  depend 
only  on  transverse  coordinates  and  fulfill  the  boundary  conditions  on  the  guide  periphery.  The 
function  expansion  is  done  in  two  dimensions  for  each  slice  separately  while  the  variations  in  the 
third  spatial  dimension  and  time  are  handled  using  the  finite  difference  approximations.  Suppose 
the  structure  was  divided  into  K  slices  in  the  r  direction  and  the  slices  are  uniformly  spaced  by 
the  distance  Ad.  Using  the  finite  difference  approximation  of  derivatives  in  the  z  directions  the 
Maxwell’s  equations  can  be  written  in  the  following  operator  form 

■^fk  =  +  L{,(st+i -yt) 

^9k  =  L|,A-1-L‘.(/*-A-i)  (1) 
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where  fk  and  gu  are  vector  functions  representing  the  electric  and  magnetic  field  at  the  ib-th  slice 
2Uld 


Lj,  =  - rV,  X  (■)  L*,  =  — — i  X  (•) 

<o<;(*,y)  foe*(x,g)Sd 

*  -1  -  -1  . 


I'Ji  = - *7 - (•)  = - *7 - TT:7i  x  (■) 

In  the  above  equations  by  i  we  have  denoted  a  unit  vector  in  the  z  direction. 

The  fields  are  now  expanded  separately  on  each  slice  according  to 

/*(<)  =  5^au(0/u(*.y)  k=l...K 


(2) 


(3) 


where  are  the  expansion  coefficients  (time  dependent)  and  /u(*,y),ffi»(x,y)  are  the  basis 

functions.  The  next  step  is  to  introduce  the  expansion  (3)  into  (1)  and  replace  time  derivatives 
by  finite  difference  formulas.  This  results  in  equations  in  which  the  only  unknowns,  for  a  fixed 
time  instant,  are  the  expansion  coefficients  at  all  slices  of  the  structure.  To  evaluate  the  expansion 
coefficients  we  take  the  inner  product  of  equations  valid  for  a  given  slice  with  expansion  functions 
and  use  the  orthogonality  property  of  the  expansion  functions.  As  a  result  we  arrive  at  [4]; 

a"  = 

jn+l/2  _  (4) 

where  Af  is  the  time  step,  a  and  6  are  column  vectors  containing  expansion  coefficients  for  all 
slices  and  superscript  n  denotes  the  time  step.  The  matrices  4  and  g  contain  the  inner  products 
and  have  the  following  str-’oture 

4  =  qdiag[4:*,4;i*]  g  =  qdiag  (5) 

The  elements  of  the  submatrices  are  given 

~  <  (^l«  “  > /<»  > 

=  <(L‘.  +  L*,)/^.,,.,  >  =  >  (6) 

Equations  (4)  show  that  expansion  coefficients  for  all  slices  are  updated  at  eeu;h  time  step  as  a 
result  of  mutual  interactions  of  fields  due  to  the  inhomogeneity  introduced  by  space  dependence 
of  constitutive  parameters.  As  far  as  numerical  cost  is  concerned,  the  most  critical  point  in  the 
PEE  algorithm  is  the  csdculation  of  the  inner  products  on  inhomogeneous  slices.  However  for 
inhomogeneous  slices  a  classical  FDTD  algorithm  can  be  used.  Combining  these  two  time  domain 
techniques  we  create  a  hybrid  method  in  which  different  algorithms  are  used  in  different  parts  of 
the  computational  space.  The  FDTD  is  used  in  the  regions  in  which  a  fine  resolution  of  field  is 
necessary  (eg.  near  edges,  media  interfaces)  and  the  PEE  is  applied  in  the  homogeneous  subregions. 
This  hybrid  approach  results  in  savings  in  numerical  effort  and  computer  memory.  This  is  because 
the  PEE  is  extremely  efficient  for  homogeneous  slices  as  matrices  *nd  are 

diagonal  so  that  the  computations  are  fast  especially  when  the  expansion~inrnction8  are  chosen  in 
such  a  way  that  they  constitute  a  set  of  eigenfunctions  of  the  Laplace  operator  defined  on  2D  region 
forming  a  slice.  In  that  case  each  expansion  function  satisfies  the  boundary  condition  and  field 
equations  globally  over  entire  slice.  As  a  result  very  few  expansion  terms  are  needed  to  accurately 
describe  field  at  each  slice.  The  FDTD  and  PEE  algorithms  are  interfaced  at  a  common  slice.  The 
transition  form  the  FDTD  to  PEE  is  done  in  the  following  way.  Given  a  field  distribution  at  the 
z  =  Zk,  provided  by  the  FDTD  part  of  the  algorithm,  the  expansion  coefficients  at  this  slice  are 
found  by  taking  the  inner  product  with  each  basis  functions  of  the  PEE.  To  switch  from  PEE  to 
FDTD  the  series  (3)  are  calculated  lA  the  interface  plane  at  the  points  required  by  FDTD. 

NUMERICAL  EXAMPLE 
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In  order  to  verify  the  hybrid  algorithm  the  cutoff  frequency  of  the  EHn  of  a  rectangular  guide 
loaded  with  a  dielectric  slab  shown  in  fig.2  was  computed  and  compared  with  the  results  obtained 
with  a  classical  FDTD  technique.  For  both  algorithms  the  identical  excitation,  space  discretization 
{Ad=  .5mm),  time  step,  and  number  of  samples  were  assumed.  The  results  are  given  in  Table  1. 
For  the  FDTD  algorithm  alone  the  CPU  time  for  the  assumed  discretization  mesh  was  53s.  For  the 
hybrid  algorithm,  the  slab  region  was  treated  with  the  FDTD  and  the  lateral  homogeneous  regions 
were  calculated  with  PEE.  The  CPU  time  depends  on  the  localization  of  the  interface  plane  and 
the  number  of  expansion  functions  used  in  the  PEE  part  of  the  algorithm.  When  the  interface  is 
at  xi  =  7.5  and  xj  =  12.5mm,  ie.  when  the  FDTD  mesh  is  terminated  only  one  slice  away  from 
the  inhomogeneous  region,  the  CPU  time  of  the  hybrid  algorithm  varies  from  18s  (for  1  expansion 
function)  to  27s  (for  5  expansion  functions),  with  most  of  the  time  (148)  consumed  by  the  FDTD 
computations.  The  error  introduced  by  low  number  of  expansion  function  is  the  Izirgest  if  only  one 
term  is  used  but  for  this  structure  is  less  then  0.1%  compared  the  result  obtained  from  pure  FDTD 
calculations,  interface  plane.  The  error  decreases  as  the  interface  plane  is  moved  away  from  the 
inhomogeneity.  This  is  due  to  the  fact  that  higher  order  terms  in  the  field  expansion  correspond  to 
higher  order  waves  traveling  in  the  lateral  direction.  If  the  interface  planes  are  located  at  xi  =  5 
and  X3  =  15mm,  only  one  term  in  the  PEE  part  is  sufficient  to  obtain  exactly  the  same  results 
as  with  purely  FDTD  technique.  The  CPU  time  for  this  case  is  31s  of  which  28s  is  spent  in  the 
FDTD  part. 


CONCLUSIONS 

A  new  hybrid  PEE^FDTD  algorithm  for  the  time  domain  analysis  of  electromagnetic  waves  in 
shielded  structures  was  introduced.  The  obtained  results  indicate  that  it  is  possible  to  obtain  the 
acceleration  of  time  domain  calculation  by  using  the  PEE  algorithm  in  homogeneous  parts  of  the 
structure.  The  results  presented  in  this  letter  indicate  that  it  is  possible  to  obtain  improvement  in 
the  speed  of  time  domain  computation  of  3D  and  2D  shielded  structures  in  which  the  homogeneous 
regions  are  predominant  such  as  microstrip  lines,  coplanar  guides,  and  discontinuities  in  planar 
guides. 
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Table  1:  Comparison  of  the  results  and  CPU  times  (DELL  486/66)  for  the  cutoff  frequencies  of 
EHn  mode  in  the  20  by  6  mm  rectangular  guide  loaded  with  a  dielectric  slab  Cr  =  2.5,  w  =  4mm. 
h  =  4mm  using  a  hybrid  algorithm  with  number  of  expansion  functions  N  and  the  localization  of 
the  interface  planes  as  pwameters 


PEE-FDTD 

Error 

rel.  to  FDTD 

■Sill 

njggilil 

CPU 

combined 

Speed  np  rel. 
to  FDTD  (538) 

Interface  of  algorithms  at  z  =  7.5,  z  =  12.5  mm 

m 

20.2075  GHz 

+  0.1% 

148 

48 

18s 

2.9 

ai 

20.1975  GHz 

148 

98 

23s 

2.3 

ai 

20.1975  GHz 

14s 

13s 

278 

1.96 

Interface  of  algorithms  at  z  =  5,  z  =  15  mm 

1  11  20.1875  GHz 

0% 

288 

3s 

31s 

1.7 

0^ 

28s 

6s 

34s 

1.55 

Figure  1;  A  structure  discretized  along  one  coordinate 


Figure  2:  Geometry  of  the  guide  used  in  the  numerical  test  showing  FDTD  and  PEE  meshes  (Slab 
centered  with  respect  to  x=a/2,  Cr  =  2.5,  dimensions:  h  =  4mm,  w  =4mm,  a=20mm,  b=6mm) 
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ABSTRACT 

This  contribution  presents  the  derivation  of  the  stability  condition  for  various  types  of  time 
domain  algorithms  used  in  the  solution  of  linear  hyperbolic  differential  equations  which  arise  in 
the  investigation  of  transient  electromagnetic  fields.  The  stability  condition  of  the  algorithm  is 
derived  by  investigating  the  properties  of  operators  in  suitably  defined  Hilbert  spaces.  Compared 
to  the  classical  von  Neumann  stability  analysis,  the  functional  analysis  approach  gives  more 
general  results  which  can  be  easily  applied  to  some  recent  and  possible  future  time  domain 
schemes. 


INTRODUCTION 

Explicit  algorithms  for  the  solution  of  initial  value  problems  have  recently  received  much 
attention  among  researchers  involved  in  the  numerical  analysis  of  electromagnetic  fields.  Two 
methods  belonging  to  this  class,  known  as  finite  difference-time  domain  (FDTD)  and  transmis¬ 
sion  line  matrix  (TLM)  algorithms  have  intensively  been  developed  in  the  last  decade.  Their 
salient  feature  is  that  electromagnetic  field  is  analyzed  in  the  time  domain  and  the  samples  of 
relevant  physical  quantities  at  nodes  located  at  the  discrete  points  in  space  are  used  to  represent 
a  physical  continuum.  These  two  methods  are  constantly  being  improved.  The  improvements 
include  the  application  of  graded  meshes  or  non  orthogonal  cells,  application  of  local  approxi¬ 
mations  or  extension  of  the  basic  algorithms  to  the  new  class  of  materials  such  as  ferrites  or 
dispersive  media.  Also  new  concepts  of  space  representation  of  fields  have  been  introduced. 
Recognizing  the  progress  achieved  in  the  recent  years  in  the  time  domain  analysis  of  electro¬ 
magnetic  fields,  it  should  be  noted,  that  the  explicit  edgorithms  underlaying  these  methods  are 
not  unconditionally  stable  and  the  improvements  introduced  to  rdgorithms  affect  their  stability. 
Consequently  there  is  a  need  to  investigate  the  stability  criteria  for  new  schemes  [6].  In  this 
contribution  we  shall  present  how  stability  of  different  algorithms  can  be  investigated  using 
functional  analysis. 

STABILITY  ANALYSIS  OF  EXPLICIT  TIME  DOMAIN  ALGORITHMS 
Let  us  consider  a  hyperbolic  differential  equation 

^/  +  I./  =  0  (1) 

where  I>  is  a  linear  differential  operator.  The  hyperbolic  equation  of  this  type,  supplemented  by 
conditions  at  <  =  0  can  be  solved  for  t  >  0  using  a  classical  finite  difference  explicit  algorithm 
[2].  It  is  a  known  fact  that  the  explicit  algorithms  are  conditionally  stable.  The  approach 
most  frequently  used  to  derive  the  stability  condition  is  known  as  the  von  Neumann  stability 
analysis  [2,  5].  This  analysis  involves  local  expansion  of  unknown  functions  into  Fourier  series 
and  assumes  the  finite  difference  representation  of  the  operator.  If  L  is  a  Laplacian,  the  von 
Neumann  approach  leads  the  formula  known  as  the  Courant  -  Friedrich  -  Levy  condition.  More 
general  results  can  however  be  derived  using  an  alternative  approach  based  on  the  functional 
analysis. 
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To  investigate  the  stability  of  explicit  algorithms  for  the  hyperbolic  equations  it  is  useful  to 
present  a  problem  in  a  canonical  form: 

AhR^f  +  \f  =  Q  (2) 

The  time  marching  algorithm  for  the  above  problem  is  stable  if  the  following  conditions  are 
fulfilled  [1]  (p.386): 


A  =  A*  >  0,  R  =  R*  >  0 

R-^  >0 
4 


(3) 

(4) 


In  other  words,  for  the  time  marching  algorithm  to  be  stable  it  is  sufficient  that  both  operators 
A  and  R  are  self  adjoint  and  positive  and  additionally  the  operator  R  —  0.25A  is  nonnegative. 
The  canonical  form  (2)  is  obtained  from  (1)  by  simply  writing  it  as 


(5) 


Where  I  is  the  identity  operator. 

Comparing  (5)  with  (2)  we  get  R  =  I/A^t  and  A  =  L.  It  can  readily  be  verified  that  operator 
L  is  symmetric  and  positive.  It  suffices  to  verify  the  condition  (4).  This  condition  is  fulfilled 
when 


\2-\\>m 

'A2<"-  4 


(6) 


or 


At< 


y/m 


(7) 


Thus  the  maximal  time  step  in  explicit  time  domain  algorithms  considered  here  depends  on  the 
norm  of  the  operator  L. 

For  a  self-adjoint  bounded  operator  L  defined  in  the  Hilbert  space  7i  the  norm  is  defined  as 

[3] 

||L||=:  sup  I  <  Ly,y  >  I  =  |A„„(  (8) 

yew.llvINi 


where  A^ar  is  the  largest  eigenvalue  of  L. 


STABILITY  ANALYSIS  FOR  ONE  DIMENSIONAL  PROBLEM 


One  important  conclusion  which  follows  from  the  functional  analysis  approach  is  that  the 
stability  condition  depends  on  how  the  unknown  functions  are  represented.  This  is  because  the 
norm  of  the  operator  depends  on  the  space  it  acts  in.  When  solving  a  particular  problem  we 
choose  the  way  the  functions  are  represented  and  the  criteria  to  measure  the  accuracy  of  our 
solution.  This  choice  is  equivalent  to  the  choice  of  a  functional  space  and  affects  the  norm  of 
the  operator  and  thus  the  stability  condition.  To  illustrate  this  problem  in  more  detail  let  us 
consider  a  one  dimensional  problem 


/(*,  fo)  =  /o(*),  fix  =  0)  =  /(x  =  /)  =  0  (10) 

where  6(x)  >  0  is  a  time  independent  continuous  function  of  x, 

One  possible  way  of  solving  the  above  problem  is  to  use  a  classical  finite  difference  approach 
but  let  us  find  the  solution  by  means  of  the  method  of  moments.  Let  D  denote  the  domain 
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of  operator  L  and  assume  that  it  allows  only  square  integrable  functions  satisfying  Dirichlet 
conditions  at  both  ends  of  the  interval.  By  equipping  the  domain  V  with  an  inner  product 


<u,v  >= 


we  specify  it  in  terms  of  the  Hilbert  space. 
It  can  easily  be  verified  that  operator 


L  = 


(11) 


(12) 


is  positive  and  self  s^ljoint.  However  if  we  would  like  to  calculate  its  norm  in  this  space  we  note 
that  the  operator  L  is  unbounded  and  consequently  its  norm  is  infinite.  Its  norm  becomes  finite 
however,  if  the  operator  is  allowed  to  act  in  a  finite  dimensional  space.  This  is  what  happens  in 
praictice  because  we  always  look  for  approximate  solution  to  the  problem  using  a  finite  number 
of  elements  to  represent  a  function.  Let  us  now  expand  the  function  f{x)  into  series  of  basis 
functions 

and  use  the  inner  product  (11)  to  find  the  expansion  coefficient  at  any  instance  of  time. 

The  finite  set  of  basis  function  defines  the  approximate  finite  dimensional  subspace  of  original 
domain.  Consider  the  following  truncated  set  of  basis  functions 


i  <  Nm 


(14) 


The  basis  functions  (14)  span  a  finite  dimensional  space  CV  in  which  the  approximate 
solution  is  sought  for.  Now  it  is  easy  to  find  the  upper  bound  of  the  operator. 


||L||  =  ||6(x)^(-)||  <  IIW~(  )||  =  ||L„ 


(15) 


Where  6mo*  is  the  maximal  absolute  value  of  b{x)  over  the  interval  <  0,/  >.  The  eigenvalues 
Ai  of  operator  Lm  are  given  by 


and  consequently  the  norm  of  L 


This  leads  to  the  condition 


~  b  /2 

(16) 

-2 

l|L|l<^ 

(17) 

At  < 

-  rrNu 

(18) 

If  we  chose  an  alternative  way  and  represent  the  function  and  the  operator  in  a  finite 
difference  sense  by  specifying  their  values  at  discrete  points  the  norm  will  be  changed.  If  the 
discretization  points  are  equidistant  and  the  spacing  is  Ad  then  [1] 


l|L||< 


{Ad)H„, 

yielding  the  Courant  -  FViedrich  -  Levy  condition; 

At< 


-  Ad 

APPLICATION  TO  ELECTROMAGNETICS 


(19) 


(20) 
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The  stability  analysis  described  above  can  be  used  in  electromagnetic  problems.  Here  the  ope¬ 
rator  L  can  be  specified  as; 

L= - ^ - -Vx  Vx(.)  (21) 

(other  definitions  are  also  possible). 

As  an  example  let  us  consider  a  cube  O  with  the  dimensions  /  x  /  x  /.  In  this  region  we  seek 
an  approximate  solution  to  the  hyperbolic  equation  with  an  operator  defined  by  (21)  given  a 
finite  number  of  expansion  functions  in  the  form  of  normalized  products  of  sines  and  cosines 

sin^^  or  cos^^  hk<NM  (22) 

Let  X>  denote  the  domain  of  operator  L.  The  basis  functions  (22)  span  a  finite  dimensional  space 
WiVv  C  T>.  We  calculate  the  upper  bound  of  the  norm  of  operator  ||L||.  Note  that  j|L||  <  ULmll. 
Where 

Lm  =  (eo/io^mtn/^min)  X  V  X  (•)  (23) 

and 

fmin  =  inffr(a:,y,i),  =  inf/ir(x,y,2)  x,y,zeii  (24) 

Using  the  same  procedure  as  for  one  dimensional  case  we  find  the  norm  of  operator  L  in 

l|I-||<||L„|l<rL,-f—  (25) 

where  Vmax  =  (foPofminA'min)”^^^  is  the  maximum  velocity  for  a  plane  wave  in  the  structure. 
Using  the  above  estimation  we  get  the  following  stability  condition 


< 


2/ 

t'mo* 


(26) 


If  the  region  Q  is  a  rectangular  prism  with  the  dimensions  a  x  6  x  (  and  the  upper  bound  for 
i,  I,  m  in  the  trigonometric  expansion  functions  is  Km,Lu,Nu  then  the  condition  (26)  becomes 

At  < - 7  -  ^  (27) 

Discretizing  the  space  D  in  the  z  direction  with  the  step  Ad  and  using  Ku  x  Lm  sine  and 
cosine  biisis  functions  to  in  the  x  and  y  directions,  the  stability  condition  derived  using  (19)  is 

At  < - ,  (28) 

For  the  discretization  of  all  three  coordinates  with  steps  Az,  Ay,  Ad  we  shall  get  the  well 
known  Courant  condition  [5] 


At  < - ^ .  — .  (29) 

W^(^)»  +  (^)^  +  (^)* 

At  this  point  it  is  interesting  to  observe  that  the  derivation  described  above  provides  stability 
criteria  for  a  few  recently  published  time  domain  algorithms.  For  instance  conditions  (27)  and 
(28)  are  the  stability  criteria  for  the  Total  Eigenfunction  Expansion  and  Partial  Eigenfunction 
Equations  schemes  derived  in  [8]  (for  sine  and  cosine  basis  functions).  In  a  compact  2-D/FDTD 
algorithm  described  in  [4]  and  investigated  subsequently  by  Cangellaris  [6],  the  functions  are 
represented  by  samples  at  uniformly  discretized  cartesian  coordinates  z,  y  and  the  variation  in 
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the  z  direction  given  in  the  form  exp{—jJz).  For  this  case  the  functional  analysis  approach 
gives 

At  < - 

^max 

This  condition  is  identical  as  the  one  given  in  [4]  and  [6].  Also  the  stability  of  a  hybrid  spec- 
tral/FDTD  method  introduced  recently  by  Cangellaris  et  al.  [7]  follows  from  condition  ,7 f. 

CONCLUSIONS 

The  application  of  the  functional  analysis  to  the  investigation  of  the  stability  of  time  domain 
algorithms  has  been  presented.  It  was  shown  that  the  method  can  easily  be  applied  to  the 
investigation  of  the  properties  of  novel  time  domain  schemes  for  Maxwell's  equations. 
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